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A COMPUTER PROGRAM FOR CALCULATING THE PERFECT GAS 
INVISCID FLOW FIELD ABOUT BLUNT AXISYMMETRIC BODIES 


AT AN ANGLE OF ATTACK OF 0° 

By Ernest V. Zoby and Randolph A. Graves, Jr. 

Langley Research Center 

SUMMARY 

A method for the rapid calculation of the inviscid shock layer about blunt axi- 
symmetric bodies at an angle of attack of 0° has been developed. The procedure is 
of an inverse nature, that is, a shock wave is assumed and calculations proceed along 
rays normal to the shock. The solution is iterated until the given body is computed. 

The flow field solution procedure is programed at the Langley Research Center for the 
Control Data 6600 computer and requires a core storage of 50000g. The geometries 
specified in the program are spheres, ellipsoids, paraboloids, and hyperboloids which 
may have conical afterbodies. 

The normal momentum equation is replaced with an approximate algebraic 
expression. This simplification significantly reduces machine computation time. 
Comparisons of the present results with shock shapes and surface pressure distribu- 
tions obtained by the more exact methods indicate that the program provides reasonably 
accurate results for smooth bodies in axisymmetric flow. However, further research 
is required to establish the proper approximate form of the normal momentum equation 
for the two-dimensional case. 


INTRODUCTION 

In the past few years, a large effort has been devoted to developing methods for 
solving the inviscid flow field over blunt bodies traveling at supersonic and hypersonic 
speeds. Early investigations are documented by Hayes and Probstein (ref. 1), and some 
of the more recent work is discussed by Grose (ref. 2). Except for the time asymptotic 
technique, the methods are primarily restricted to the subsonic and slightly supersonic 
part of the flow field with the downstream supersonic flow being computed by the method 
of characteristics. A simple approximate but accurate method was developed by Maslen 
(ref. 3) and was shown to be applicable throughout the subsonic and supersonic regions 
of the flow field. This method was extended to include nonequilibrium chemistry in 
reference 2 where it was shown to provide significant reductions in machine storage 
and computing time over the more exact methods. However, the method as developed 
in reference 3 is inverse, that is, the shock shape is specified and a resultant flow 


field and body are computed. This situation is undesirable since the body shape is 
usually specified and the flow field and shock shape are to be determined. In addition, 
the method neglects the normal velocity component and its gradients in the normal 
momentum equation; these factors account for recompression across the shock layer 
in the stagnation region. Jackson (ref. 4) devised an iterative scheme to solve the 
direct problem based on the analysis of reference 3. 

In this paper an approximate method has been developed for solving the direct 
problem of a perfect gas inviscid flow field over blunt bodies at an angle of attack of 0°. 
The solution incorporates the normal velocity component and the effect of its gradients 
by replacing the normal momentum equation with an analytic expression developed by 
Maslen (ref. 5). The present solution uses the basic iterative technique of reference 4 
which scales the shock to the specified body in the subsonic and low supersonic region 
of the flow field. In the downstream supersonic section of the flow field, the shock 
shape is computed by a "marching" procedure for successive points. The method 
provides a rapid technique for computing the perfect gas inviscid axisymmetric flow 
field about blunt bodies such as spheres, paraboloids, ellipsoids, and hyperboloids 
with an option for conical afterbodies. 

This paper presents a discussion of the computational method and results and 
will serve as a user's manual for the program. A comparison of the results with the 
results of more exact methods indicates that this method is reasonably accurate. The 
method is programed at Langley Research Center for the Control Data 6600 computer 
and requires a core storage of 50000g. A detailed discussion of the program is given 
in appendix A, and the program listing with a sample case is given in appendix B. 

SYMBOLS 


Symbol 


Program name 

a 

nondimens ional major axis (parallel to flow) of 
ellipse, a*/« 

ABA 

b 

nondimensional minor axis (perpendicular to flow) 
of ellipse, b*/t 

ABB 

d 

nondimensional diameter, d*/ fi 


Ht 

* 

H t 

nondimensional total enthalpy, — - — 

ur/2 

HSO 

h 

h* 

nondimensional static enthalpy, .2 / 

U*/2 



2 



Symbol 


Program name 


J denotes step through shock 


L 

body type 

L 

8 

characteristic length (all dimensions 



nondimens ionalized with this term) 


M s 

shock Mach number 

MS 


free -stream Mach number 

MINF 

N it 

iteration number 

IT 

N t 

number of points across shock layer 

NTl 


not including shock point 


N x 

number of rays in subsonic region 

NX 

n 

flow type (0 = Two dimensional; 

N 


1 = Axisymmetric) 


p(x, y) 

nondimens ional pressure through shock 
layer, 

(p*U*2 )oo 


P(^) 

nondimens ional pressure through shock 

, P*U) 

layer, — 

(p^* 2 ^ * 

PI 

Pb 

Ph 

nondimens ional body pressure, bL_ — 

(pnj* 2 )^ 

PI(I, NTl) 

Ps 

Ps 

nondimens ional shock pressure, — „ — 

( P*!!* 2 )^ 

PS 

Pst 

nondimens ional body stagnation-point 

PST 


Pst 

pressure, — 


* 

Poo 

free -stream pressure 


AR 

increment in radius of curvature 

DR 
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Symbol 


Program name 


nondimens ional body radius of curvature at 
stagnation point, R^A 

R s nondimensional shock radius of curvature, RgA 

R s t nondimensional stagnation-point shock radius 

* / 

of curvature, R s tA 

r nondimensional cross-sectional radius, r*A 

r^ nondimensional cross-sectional body radius, r^yv 

r c analytically computed nondimensional cross- 

sectional body radius, r *A 

r^ nondimensional maxitnum cross-sectional 

shock radius, r^/fi 

r n nondimensional nose radius, r */k 

r g nondimensional cross-sectional shock 

radius, rgA 

S nondimensional entropy 

jjg 

free -stream velocity 

u nondimensional shock layer velocity in 

x-direction, u*/U^ 

u s nondimensional velocity behind shock in 

x-direction, u|/U^, 

v nondimensional shock layer velocity in 

y-direction, v*/U<J, 

vs nondimensional velocity behind shock in 

y-direction, vg/U*, 


RKB 

RKSX 

RKSl 

RI 

RI(I,NT1) 

RW4C 

REND 

RS 

UI 

US 

VI 

vs 
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Symbol Program name 


X 

nondimensional coordinate measured along 
shock wave from axis of symmetry, x*/ e 

XS 

x t 

nondimensional total arc length to limiting 
Mach number, xj'/e 

XTC 

y 

nondimensional coordinate normal to shock, y*/t 

YI 


nondimensional value of y from shock to body 

YI(I, NTl) 

z 

nondimensional distance measured along axis 
of symmetry from shock wave (see fig. 1) 

ZI 

zb 

nondimensional coordinate measured along 
axis of symmetry, z£A 


z s 

nondimensional value of z for shock point, 
z|/« 

ZS 

7 

ratio of specific heats 

GAM 

<5 

nondimensional shock stagnation-point 
standoff distance, 6 *A 

ZST 

6 

convergence criterion, 1 X 10"^ 

EM 

e 

error 

E4(IT) 

AO 

increment in shock angle 


e b 

body angle 

THB 

d s 

shock angle 

THS 

p 

nondimensional local density, p*/p* x 

RHI 


nondimensional shock density, Ps/pto 

RHOS 

* 

P oo 

free -stream density 
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i Symbol 


Program name 

0 st 

second derivative of stagnation -point shock 
radius of curvature 

SDERRX 


nondimensional local stream function through 
shock layer, ^*/p£ U£«2 

SY 

^s 

1 

nondimensional shock stream function, 
4 >*/ p * U*^ 2 

o/ oo CO 

SFS 


An asterisk * denotes a dimensional quantity. An n as a subscript denotes Nr. 

RESULTS AND DISCUSSION 


The inviscid flow field over blunt bodies is described by the governing equations 
for an inviscid, nonconducting gas. The equations, for this analysis, are written in 
an orthogonal curvilinear coordinate system fixed at the axis of symmetry of the shock. 
(See fig. 1. ) Assumptions that the shock layer is thin and that the flow field is axisym- 
metric are made so that the resulting equations are 

Continuity: 


/ b (v n p u)\ 

V ax /y 


/9(r n pv)\ 

V )x 


0 


x-momentum: 


u 




uv 1^ /3p\ 

R s “ p \ax/y 


y-momentum: 


u 





Entropy: 




0 


( 1 ) 


( 2 ) 


(3) 


(4) 
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Analysis 


By following the analysis of reference 3, a Von Mises transformation is introduced 
so that the independent variables are the distance along the shock x and the stream 
function ^ . The stream function is defined by 


9i p 

ay 


r n pu 


(5) 


9ij/ 

9 x 


r n pv 


( 6 ) 


If these transformations are applied to the governing equations, the continuity equation 
is automatically satisfied; the y-momentum equation becomes 


r n9p = _u + /3v\ 
3* R s \3X/* 


(7) 


and the entropy equation is 



(8) 


The two momentum equations and the condition that the flow is isoenergetic implies 
that 


H t = h + u 2 + v 2 

In addition to equations (7), (8), and (9), the equation of state is 
h = h(S,p) 


or 


P = p(S,p) 
For a perfect gas 


h 


27 p 
7- 1 P 


(9) 


( 10 ) 


(ID 


( 12 ) 
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and for the condition of isentropic flow along a streamline (eq. (8)) 


P 



(13) 


Note that all the previous parameters were presented in a nondimensional form. Dis- 
tances are divided by a characteristic length « ; density and velocities by their respec- 
tive free-stream values; pressure by ^>*U*^) ; the enthalpy by and the stream 

function by (p*!!*)^ e n+1 . 


In the analysis of reference 3, Maslen assumed that in equation (7) is 

small and that the streamlines are nearly parallel to the shock, that is, v« u. Thus, 
equation (7) can then be approximated as follows with u and r replaced by u s and 
r s , respectively, 


9^ r n R 
x s 

Integrating from the shock toward the body yields 

P (*) = P s + (* - ^ 

r2 R 0 


(14) 


(15) 


This analytic equation allows an independent evaluation of the pressure throughout the 
shock layer and removes many of the usual mathematical problems associated with the 
subsonic and supersonic regions. 


With equation (15) and an equation for the shock shape, a simple procedure can be 
set up for computing the flow field. The calculations proceed along rays normal to the 
shock until the body = 0) is reached. Maslen (ref. 3) concluded that the accuracy of 
such a technique is adequate for general purposes. However, as previously noted, the 
analysis neglects the normal velocity term in equation (9) and the term in the 

j term contributes to the recompression in 

the stagnation region. In this investigation, preliminary calculations verified that the 
body pressure distributions resulting from equation (15) were approximately 10 to 15 per- 
cent lower in the stagnation region than results obtained from more exact methods. 


normal momentum equation (7). The (|^ 


In a more recent investigation, Maslen (ref. 5) considered the term and 

by following an analysis similar to reference 3 developed an approximate analytic 
expression for the pressure. This expression can be written for axisymmetric flow as 
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pU) = P s + 


u c 


r n R 

L s 


(l/< - pj + 


v s tan 0, 
In 


cos e\\p -p. 


Rc 


2 ^< 


(16) 


Equation (16) is the pressure equation programed in the present computer code. 

With the effects of the normal velocity gradients included in equation (16), an 
approximation to the normal velocity component is necessary. The constant density 
solution for a sphere (ref. 1) was used as a model, and an approximate expression was 
developed for this investigation as 

v(^) = v s — (17) 

Preliminary calculations indicated that equation (17) was satisfactory as an initial guess 
for the normal velocity profile, but better agreement was obtained with existing results 
if the normal velocity component is reevaluated during each iteration by 


v(*) 



(18) 


This equation was determined after several trial calculations and comparisons of the 
data. The rationale for equation (18) is simply the improved agreement of results with 
published data. 

The y-coordinate measured normal to the shockwave is computed with equa- 
tion (5) written as 


r n a_y _ _ 1 

dt p pU 

and from the flow field geometry (fig. 1) 


r = r g - y cos 0 g 


(19) 


(20) 


Substituting equation (20) into equation (19) and solving for y gives 


y 


COS 0 C 


/ O. r»nc ft 




( 21 ) 
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for axisymmetric bodies and 

y=^ s A dlp (22) 

for two-dimensional bodies. In addition, the distance measured along the axis of 
symmetry from the shock wave is 

z = z g + y sin d s (23) 

Although the use of equations (16), (17) or (18), and (21) with an equation for the 
shock shape may produce better agreement with published results than the method in 
reference 3, the method is still inverse, that is, the computed flow field corresponds 
to a given shock. The problem of greater engineering interest is the direct problem, 
that is, computing the shock shape and flow field for a specified body. 

Problem Description 

The basic steps of an interactive procedure presented by Jackson (ref. 4) are 
incorporated in this investigation along with the previously discussed flow equations 
for developing a direct solution of the perfect gas flow field over blunt axisymmetric 
bodies at an angle of attack of 0°. Note that although the use of the analytic pressure 
equation (eq. (16)) removes many of the usual mathematical problems associated with 
the subsonic and supersonic regions, the solution of the direct problem will require 
separate iterative schemes for the two regions. 

In the supersonic region, changes in body shape at a point cannot affect the shock 
shape or flow field at upstream points; therefore, a step-by-step marching procedure is 
used. The shock shape and flow field are computed at successive downstream points. 
However, disturbances, for example, changes in the body shape, in the subsonic region 
influence the entire region, and this fact requires an iterative procedure for describing 
the entire subsonic region at once. As discussed by Jackson (ref. 4), the need for 
computing the entire subsonic— low-supersonic region simultaneously can be illustrated 
by two flat -faced cylindrical bodies with different body radii. The shock standoff 
distance and curvature at the stagnation point are different for the two bodies since the 
scale of the respective subsonic regions is different. Consequently, the shock shapes 
cannot be calculated by only considering the body shape near the axis of symmetry and 
computing ahead. 

A summary of the iterative techniques is given in the following paragraphs. A 
more detailed discussion of the iterative techniques is included in the program descrip- 
tion in appendix A. 
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For the subsonic region, the shock shape is represented by a Taylor series in 
terms of the arc length as 

R s (x) = R st +a st (24) 

The shock shape is computed from the axis of symmetry to a point x^. on the shock 
defined by a shock Mach number of approximately 1.3. For the remainder of this 
paper, this region is referred to as the subsonic region. The conditions immediately 
downstream of the shock are determined from the Rankine-Hugoniot relations and the 
equation 


dd 


s 


dx 

r^j 


(25) 


Along rays perpendicular to the shock, the flow properties are calculated at local 
values of the stream function by equations (16), (17) or (18), (13), and (9). The calcu- 
lations continue until the body ('/'= 0) is reached. The geometry of each grid point is 
calculated by equations (21) or (22), (20), and (23). The computed body will probably 
not agree with the specified body in size or shape. The shock shape and flow field are 
then scaled to give approximate agreement with the specified body. The difference in 
the shape of the specified and computed body is defined by an error e . This error is 
the difference between the radius of the calculated body r c and the radius of the 
specified body at the last computed body location in the subsonic region. Depending on 
the magnitude and sign of the error term e , the a g ^ value in equation (24) is changed 
and a new shock shape, a new calculated body, and a new error result. For consistency 
on successive iterations, the position of xj- is computed so that it corresponds to the 
same shock angle and Mach number as on the first iteration. The process is repeated 
until the error is less than some convergence criterion (e = 10"4 for this program). 

For this region the specified bodies are given in the program by analytic equations. 
A sphere, an ellipsoid, a paraboloid, or a hyperboloid may be specified. 

If the calculation of the supersonic region is needed, a step-by-step marching- 
ahead procedure is used to compute the shock shape and flow field. The last converged 
values of the shock radius of curvature and shock angle in the subsonic region are 
incremented. The arc length is computed and held constant for the following iterations. 
By holding the arc length rather than the shock angle constant for iterations in the 
supersonic region, the shock inflection point resulting from the overexpansion of the 
flow and the subsequent recompression were found to be easier to compute. The flow 
properties and the geometry of the grid points are then computed with the same equations 
as for the subsonic region. The calculated body is compared with the specified body, and 
an error is computed. A new increment in the shock radius of curvature is determined 
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and the process is repeated until satisfactory convergence is obtained. The procedure 
is used to extend the solutions to any downstream distance r^. As the calculations 
proceed downstream of the shock inflection point, an option is included for computing 
sharp cone results (ref. 3) when the shock angle is within a small value (0. 001 radian) 
of the equivalent sharp cone shock angle. This option was included since it was difficult 
to obtain satisfactory convergence when the shock is nearly straight and the shock radius 
of curvature is extremely large. In addition to the calculation of the flow field over the 
analytic body shapes discussed previously, an option is included for computing the 
supersonic flow fields over conical afterbodies. 

Comparison of Results 

A comparison of results of existing computational methods for shock shapes and 
surface pressure distributions with results of the present method is given in figures 2 
to 11. These computational methods used for comparison are based on the solution 
of the full governing inviscid flow equations. Comparisons are presented for a range 
of Mach numbers and body shapes. 

In figures 2 to 8, a comparison of existing shock shape and surface pressure 
distributions for hemispheres, paraboloids, ellipsoids, and hyperboloids (refs. 6 to 8) 
with the corresponding results of the present method is shown. The agreement between 
the present method and the published results is, in general, very good. 

In figure 9 the shock shape and pressure distribution for a hemisphere-cylinder 
at = 7. 7 (ref. 9) are compared with the present results. The present calculations 
yielded shock standoff distances over the cylindrical section of the body which are 
approximately 12 percent smaller than those obtained from reference 9. This comparison 
is similar to that shown in reference 3. The reason for this disagreement may possibly 
be explained by the assumption in the present analysis of a thin shock layer which is not 
strictly valid over the cylindrical afterbody. However, the body pressures are predicted 
very accurately. 

In figure 10, the shock shape and pressure distribution over a spherically blunted 
30° half-angle cone for M^ = 8 and y = 1.4 computed by the method of reference 10 
are compared with the present results. The agreement is very good. At a shock axial 
location of approximately 2 (see fig. 10(c)), there is an inflection in the shock shape. 

This shock inflection results from the overexpansion of the flow and subsequent recom- 
pression to cone conditions. In figure 11, a similar comparison is made on a spherically 
blunted 40° half -angle cone at M^ = 38. 4. The flow field computed by the method of 
reference 10 is in chemical equilibrium. A good match of the shock standoff distance 
was not obtained with the present method for 7 = 1.4 as would be expected. The 
shock standoff distances for the present method shown in figure 11 were obtained with 
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7 = 1. 16. This is the shock stagnation value computed by the method of reference 10. 
Pressure distributions are shown for 7 = 1.16 and 7 = 1.4. Both agree reasonably 
well with the results obtained by the method of reference 10. 

The results presented previously were for axisymmetric blunt bodies. The 
present computer code can be applied to two-dimensional bodies. However, initial 
comparisons of results from reference 11 on two-dimensional cylinders for Moo from 
4 to 20 and results from the present method indicated some problem areas. Primarily, 
it was noted that the computed pressures across the shock layer were higher than the 
body stagnation-point pressure. In developing equation (16), Maslen (ref. 5) assumed 
the dependence of the recompression term and pointed out that equation (16) was an 
approximation whose validity could only be established by comparisons with experimental 
data or exact theoretical analyses. Since the results of the present investigation indicate 
that the recompression term suggested by Maslen is applicable to axisymmetric bodies 
but not to two-dimensional bodies, further research is required to establish the proper 
form for two-dimensional bodies. 

From the results presented, it is shown that the present program provides 
reasonably accurate estimates for quantities such as standoff distances, shock shapes, 
and surface pressure distributions for axisymmetric flow over smooth bodies. 

Existing data for the flow parameters through the shock layer were not as readily 
available as the shock shape and pressure distribution data. Qualitative comparisons 
with the results of reference 12 were good except for the profiles near the stagnation 
line. For example, the tangential velocity is more nonlinear than the results of 
reference 12 in this region. With the approximations in equations (16) and (17) or (18), 
such deviations, especially in the stagnation region where the effects of the normal 
velocity are predominant, are not surprising. In general, the present program provides 
a reasonably good flow field description. More accurate results could probably be 
obtained by refining the expression for the normal velocity component (eq. (18)). 

CONCLUDING REMARKS 

A computer program has been developed for solving the direct problem of a 
perfect gas inviscid flow field over blunt bodies at an angle of attack of 0°. The geome- 
tries specified in the program are spheres, ellipsoids, paraboloids, and hyperboloids 
which may have conical afterbodies. The calculation procedure is of an inverse nature, 
that is, a shock wave shape is assumed and calculations proceed along rays normal to 
the shock, but the shock shape is iterated until the specified body is computed. 

The simplifying assumption which significantly decreases machine computation 
time is that the normal momentum equation is replaced by an explicit analytic expression. 
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Although these calculations are of an approximate nature, a comparison of the present 
results with existing shock shapes and surface pressure distributions indicates that 
the program provides reasonably accurate results. 

The option for two-dimensional bodies is included in the program. However, 
attempts to apply the present approximate form of the normal momentum equation to 
two-dimensional bodies were not successful, and further research is required to 
establish the proper form for the two-dimensional case. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va. , August 30, 1973. 
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APPENDIX A 


DETAILED PROGRAM DESCRIPTION 


As an aid to a better understanding and a further description of the calculation 
procedure, flow diagrams are presented with a discussion of the processes. 

Flow Diagram for Subsonic Calculation Procedure 
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APPENDIX A - Continued 
Flow Diagram for Supersonic Calculation Procedure 


Increment last values for 
shock radius of curvature 
and shock angle 


Compute average radius 
of curvature 


e(x) = o s + ^ • J 
AX-j = Ravg 1 • 



ar 2 < ar 1 


Compute AR 


ar 2 > ar 1 


>s - THS C 



Compute \' 

? 

flow field 
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APPENDIX A - Continued 


Input 

The input for the program is very simple and is primarily restricted to subroutine 
INITIAL except for the cone half-angle THC required in subroutine SUPER. In 
INITIAL, the necessary inputs are the free-stream Mach number MINF, the ratio of 
specific heats GAM, the type of body denoted by L, the type of flow denoted by N, and 
the radius of curvature of the body RKB. The number of rays is given by N x or NX 
and the number of points on the rays is N^ or NTl. A maximum of seven and a 
minimum of four rays can be used. Similarly, a maximum of 15 and a minimum of 10 
(for consistent results) points along the rays can be used. The stagnation line is 
considered to be a ray although no calculation is made along this ray. Also, the number 
of points on a ray do not include the shock point. 

Shock Shape 

Equation (24) is used to represent the shock-wave shape in the subsonic region. 

The equation is calculated in the subroutine SADCURV where R s j-, the shock radius of 
curvature at the stagnation point, is given as RKS1 and a g ^, the second derivative of 
the radius of curvature at stagnation point, is given as SDERRX(IT). For the first 
iteration (IT), RKS1 is assumed to be equal to the body radius of curvature RKB, and 
SDERRX(l) is 0. The computed value of RKS1 from the procedure is used on further 
iterations. For the second iteration, SDERRX(2) is assigned an arbitrary small posi- 
tive or negative value as discussed in the section on iteration procedure in this paper 
and on subsequent iterations SDERRX(IT) is computed. 

With the shock radius of curvature and the arc lengths, the shock angle (THETAS) 
is computed by equation (25) with a six-point integration routine INT6. The shock angle 
and the corresponding pressure (PS), velocities (normal (VS) and tangential (US)) and 
density (RHOS) behind the oblique shock wave are computed at each shock location in sub- 
routine SONIC until the Mach number behind the shock (MS) is equal to approximately 1.3. 
From experience it was seen that a limiting Mach number near this value produced a 
smooth transition from the subsonic to supersonic iteration procedure. The angle from 
the stagnation point to the point on the shock defined by the limiting Mach number is held 
constant for all iterations so that only the arc length is altered on each successive itera- 
tion. The radial (RS) and axial (ZS) locations of each shock point are then determined 
from the flow field geometry in subroutine SUBGEM. However, for a more rapid pro- 
cedure, the properties across the shock layer are not computed at each of these shock 
locations. The arc length to the limiting Mach number is divided into NX - 1 equally 
spaced regions and NX points on the shock are computed in SUBGEM. The shock geom- 
etry and properties at each of these points are then determined from an interpolation 
scheme in subroutine SUBPROP. 
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APPENDIX A - Continued 


Flow Properties 

The simplifying feature of this program is that the normal momentum equation is 
replaced by an analytic expression from reference 5. This expression is given as 
equation (16) in this paper and is programed in subroutine PRESLOC. 

In references 3 and 4, the normal velocity component (VI) is assumed to be much 
less than the tangential component (UI) even in a stagnation region and is consequently 
neglected. An approximation to the normal velocity is given in subroutine PROPLOC. 

For a first guess, the equation is given as equation (17) in this paper and is reevaluated 
by equation (18) during each iteration of the shock shape. 

The local density is also computed in subroutine PROPLOC. The point on the 
shock where the local stream function crossed is determined and the density is then 
simply computed from isentropic relations (eq. (13)). With these flow properties, the 
tangential velocity is computed from the energy equation. 

The local flow properties are computed at points along each ray until the body (\p - 0) 
is reached. The points are defined by a local stream function which is computed from the 
stream function at the shock point and equally spaced increments of the stream function 
along each ray. 


Iteration Procedure (Subsonic) 

With the flow properties computed, the location (radial, axial, and normal) of 
each local stream function through the shock layer is computed in subroutine GIP. In 
order to relate these dimensions to the known body geometry, a scale factor is computed 
in subroutine BADCURV. The equation for a specified body is passed through the body 
points computed from the flow field analysis (subroutine GIP), and a relevant length (for 
instance, the body radius of curvature in the case of a sphere) is computed. The scale 
factor is the ratio of this relevant length to the known corresponding length for the desired 
body. All dimensions, shock and body, are then scaled by this ratio in subroutine SCALE. 
The scaled body points as noted in the subroutine are used with the appropriate analytic 
equation in subroutine SOSOD to compute the shock stagnation point standoff distance ZST. 
The value of ZST and the scaled axial length for the last body point in the subsonic region 
are then used again with the proper analytic body equation to calculate the radial com- 
ponent of this body point in subroutine R4E4. With the scaled and the analytically computed 
values of this radial location, an error value E4(IT) is also computed in this subroutine. 
This error value is the percent difference in the two radial calculations. For the second 
iteration, a small negative value of a st is assumed if E4(l) is positive, and a small 
positive value of <? s t is assumed if E4(l) is negative. For subsequent iterations, o s ^. 
is computed by a Newton -Raphson technique in subroutine SDCORR. A new arc length XTC 


18 


APPENDIX A - Continued 


to the limiting Mach number is then computed in subroutine CORRl or CORR2 depending 
on the sign of a s t . The problem is now iterated until desired convergence is obtained. 

Output (Subsonic) 

The converged subsonic results are the output of subroutine OUTPUT. The shock 
and body data are given at six equally spaced shock arc lengths XSP. Note that the 
complete stagnation- line results are not computed in this program. The flow properties 
and geometry at each grid point along a ray normal to shock point XSP are also given. 

Iteration Procedure (Supersonic) 

Although the shock and flow properties are computed in the same manner as in 
the subsonic region, the iteration procedure is different. No scaling is used and a 
"marching” technique is employed along the shock. The technique is basically as follows: 

(1) A small increment for the last converged values of the shock radius of curvature 
and shock angle are assumed. 

(2) An average radius for the computing interval is calculated and the shock angles 
are computed within this interval. 

(3) The change in arc length is computed and held constant for subsequent iterations. 

(4) The shock geometry and properties are computed. 

(5) The flow field geometry and properties are computed. 

(6) A test is made if a conical option is required, and if it is, another test is made 
to determine whether the juncture point between the nose and the afterbody is reached. 

If the juncture point is reached, the body radius point is computed by the routine. When 
the two tests are negative, the analytic body radius point is computed by the equation 
used in the subsonic region. 

(7) An error value is computed and a new increment for the shock radius of curva- 
ture is also computed. 

(8) The calculations are iterated until the desired convergence is obtained. 

(9) A new step is made along the shock and the same procedures as in steps (2) to 
(8) are used. 

(10) If a shock inflection point is computed, a test is made to determine whether 
the subsequent computed shock angles are within some small value (0. 001 radian) of 
the equivalent sharp-cone shock angle. If the shock angle is sufficiently close to this 
value, a sharp-cone solution is used for the following calculations. 

(11) Calculations continue until some defined end point on the shock REND is computed. 
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Output (Supersonic) 

The converged supersonic results are the output of subroutine OUTPUTS at each 
successive downstream point. The shock and body data are given for this region. 
Other data, for example, flow properties and geometry of grid points, may be written 
if desired. The symbols in subroutine OUTPUTS are defined as follows: 

ZW body axial location measured from axis of symmetry of shock 

RW body radial location measured perpendicular to axis of symmetry 

YW body normal location measured perpendicular from shock point 

ZS shock axial location measured from axis of symmetry of shock 

RS shock radial location measured perpendicular to axis of symmetry 

PS pressure at shock 

RHOS density at shock 

US tangential velocity at shock 

PW body pressure 

RHOW body density 

UW body tangential velocity computed normal to ray from shock 

MW Mach number at body 

THB body angle, 0 ^ , radians 
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PROGRAM LISTING AND SAMPLE OUTPUT 


4 

5 

6 

7 

8 


9 


The program listing and sample output follows: 


PROGRAM DIP GAP 

THIS PROGRAM COMPUTES THE PERFECT GAS FLOW FIELD ABOUT A BLUNT TWO 
DIMENSIONAL OR AXl SYMMETRIC BODY. THE BODY IS KNOWN AND AN INVERSE 
TECHNIQUE IS ITERATEO UNTIL OESIREC BODY IS FOUND. 

THIS MAIN ROUTINE CONTROLS THE SEQUENCE OF CALCULATIONS 

REAL MINFt MS IN2 , MS ,M SO ,MSI 

COMMON / BLK L/ GAM, M IN F, X IN IT AL , 1 1 N l T AL «R INI T AL , THETA IN , SF I N I TL , L , N 
I, SDK XI NT, RKBN, ABA, ABB, PI ,HPA,HPB 

COMMON /BLK 2/ I T , DEL X , I BEG , EM, NC , NT! «NX ,REND, NO , P ST , OR ( 61 ,F (2001 . 1 
IP 

COMMON / BLK 3/ X SP ( 7) , XT , XTC , SCP ,R KS , RK SC , RKS I , SDERRX { 8 ) ,H SO , XS1 
COMMON /BLK4/ XS(300) , RKSX (300 ) , ZS ( 3 00 ) , RS ( 300 ) , S FS 4 300 ) , Tb ETA S ( 30 
1 0 ) , US( 300 ) , PS ( 300) ,RHDS(300J , MSI 300 ) , VS ( 300) 

COMMON /BLK 5/ R I ( 7 , 15 ) , Z 1 1 7, 1 5 ) , Y I ( 7 , 1 5 > , R I C I 7 , 1 5) , L IC ( 7, 15 ) , Y IC 4 7 
1,15),RW4C,ZST,E4( ti),THB( 50) 

COMMON / B LKo/ SFC(7) , PO (7), 20 (/), PQ( 7) ,UQ( 7), TH£TAO( 7), RH0S04 7) ,MS 
iO( 7 ) , RkO 17), DEL S Y ( 7) ,RK0C(7) ,ROC( 7 ) , ZCC ( 7 ) , V J ( 7 ) 

COMMON /BLK 7/ S Y ( 7 , 15 ), P I ( 7 , i 5 ), R HU 7, 1 5 ), U I ( 7 , 15 ), MS I ( 7, 15 ), V I ( 7 , 
115) 

CALL INITIAL 
I T= 1 ♦ I BEG 
CALL SAOCURV 
CALL SONIC 
CALL SUBGEM 
CALL SUBPRCP 
CALL GIP 

IF ( IP.EQ.2) GC TO 3 

IP*2 

GO TO 2 

CALL B ADC UR V 

CALL SCALE 

CALL SOSOD 

CALL R4E4 

IF ( AB S ( E 4 ( IT) ) . LE • EM I GC TO 7 

THE SUBSONIC REGION IS COMPUTED IF THIS CRITERIA IS MEET. NCR GO 
TO EXIT OR SUPERSONIC REGION. 

IP«1 

IF ( IT.EQ.1.AND.E4(1) .GT.O.OJ ^0 TO 5 
IF (IT.GT.l) GO TO 8 
l T * I T* L 

SDERRX ( 2) =2. 

CALL CORK 1 
GU TO 1 
I T= I T+ 1 
S DERRX ( 2 ) --2 • 

CALL C0RR2 
GO TU 1 
CALL OUTPUT 
GO TO 9 
CALL SDCURR 

IF 4 SDERR X ( IT) .GT.O.O) GC TO 4 

GO TO 6 

CALL SUPER 

STOP 

END 


A 1 
A 2 
A 3 
A 4 
A 5 
A 6 
A 7 
A 8 
A 9 
A 10 
A 11 
A 12 
A 13 
A 14 
A 15 
A 16 
A 17 
A 18 
A 19 
A 20 
A 21 
A 22 
A 23 
A 24 
A 25 
A 26 
A 27 
A 28 
A 29 
A 30 
A 31 
A 32 
A 33 
A 34 
A 35 
A 36 
A 37 
A 38 
A' 39 
A 40 
A 41 
A 42 
A 43 
A 44 
A 45 
A 46 
A 47 
A 48 
A 49 
A 50 
A 51 
A 52 
A 53 
A 54 
A 55 
A 56 
A 57- 
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subroutine initial b i 

B 2 

THIS ROUTINE INITIALIZES QUANTITIES PLLS DEFINES CONSTANTS B 3 

B 4 

REAL MINF,MSIN2.MS,MS0,MSI B & 

COMMON /bLKl/ GAM,MINF,XINITAL,ZINITAL,RIMTAL,THETAIN,3FINITL,L.N B 6 

l.SDRXINT.RKdN, AE A , ABe , P 1 , PPA , HPB B 7 

COMMON /&LK2/ IT,0£LX,lBEG,EM,NC,NTl,NX,REND,ND,PST,DRl8),F(200i, I B 8 

IP B 9 

CCMMCN /BLK3/ XS P I 7 ) , XT , XTC , SC P , RkS » RK SC , RKS i , SDERRX ( B J , H SO , XSi B 10 

IP 3 1 8 11 

M INF =10 . B 12 

GAM=1.4 B 13 

XINITAL=0.0 B 14 

ZINITAL=0.0 B IS 

RINITAL 3 0.0 B 16 

THETAI N=1.57C 7963268 B 17 

SPIN ITL=0 .0 B 18 

IBEG=0 B 19 

NT 1= 15 ' 8 20 

RENU-1.5 S 21 

NX=7 B 22 

SDRxINT=0 . 8 23 

DELX 3 .015 B 24 

EM 3 . 0001 B 25 

8 26 

MINF IS THE FREE STREAM MACH NUMBER; GAM*RAT10 OF SPECIFIC FEATS B 27 

DELX 3 IN l T I A L ARC LENGTH I NCRMENT ; NT 1 =NUHdER UF POINTS THRU FLOW B 28 

FIELD FOR NUMBER OF R A YS I NX- 1 )NOR M AL TO SHOCK ;REND 3 RADIUS OF SHOCK B 29 

WHERE COMPUTATIONS STOP ; EM=CONVEPGENCE CRITERION B 30 

B 31 

THH 3 70. B 32 

T HH 3 0.01 745329252*THri B 33 

HPB 3 1 . B 34 

HPA=HPB/TAM THH) B 35 

ABA=4. B 36 

ABB 3 1 . B 37 

P1=1.0 8 38 

RKBN= 1 . B 39 

RKB*P1 B 40 

RKS 3 KK8 B 41 

B 42 

RKB IS RADIUS OF CURVATURE OF BODY. FOR SPHERICAL NOSE » RKB EQUALS B 43 

RADIOS OF SPHERE (RKBN) . FOR PAR ABOL A , RKe EUCALS DISTANCE FROM B 44 

FOCUS TO DIRECTRIX(Pl) . FOR ELLIPSE. RKB ECUALS MINOR AXIS SQUARED B 45 

RATIOEO TO MAJOR AXIS! ABB**2/ ABA) . FOR HYPERBOLA , RK B EQUALS MINOR fl 46 

ASYMPTOTE SQUARED RATIOED TU MAJOR A SY MPTOTE I HPB**2/HP A J . THH EQUAL 6 47 

S HALF AKoLE OF HY PER BOLA. NOT E ... MAJOR REFERS TO LENGTH PARALLEL B 48 

TO ONCOMING STREAM. B 49 

8 50 

L=3 8 51 

B 52 

IF L EQUALS 1. THE COMPUTED RADIUS OF CURVATURE OF SPHERE h ILL BE B 53 

SCALING PARAMETER. IF L EQUALS 2, THE COMPUTED MAJOR AXIS OF ELLIPS B 54 

E WILL BE SCALING PARAMETER. IF L EQUALS TO 3, THE COMPUTED DISTANC B 55 

E FROM FOCUS TO DIRECTRIX OF PARABOLA WILL EE SCALING PARAMETER. B 56 

if l equals 4 , the computeo major asymptote cf hyperbola will be b 57 

SCALING PARAMETER. 8 58 

B 59 

X 8 CO 

B 61 

IF N = l, THE BODY IS A XI S YMMETR I C. IF N = 0, TFE BODY IS 2-D. B 62 

8 63 

RETURN 8 64 

END 3 65- 
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SUBROUTINE SAOCURV 

THIS ROUTINE COMPUTES THE SHOCK RADIUS OF CURVATURE AS FUNCTION OF 
X. RKS1 IS THE STAGNATION POINT VALUE. RKSX IS VALUE ALONG SHOCK. 
XS=SHOCK ARC LENGTH 

COMMON / BLK 1/ GAM, MINF, XI NIT AL, ZI NIT AL ,RINI TAL ,THETAIN # SF I NITL *L ,N 
l, SDR XI NT,RK BN, ABA, ABB, PI ,HPA,HPB 

COMMON / BLK2/ I T , DELX , I BEG, EM , NC , NT 1 ,NX, REND , NO , PST , DR ( 8) , F (200) , I 
IP 

COMMON /BLK 3/ XSP ( 7) , XT , XTC, SCP ,RKS , RK SC ,RK$1 , SOERRX l 8 1 ,H SC ,XS1 
CQMMCN /BLK4/ XS (300 ), RKSX (300) ,ZS (300) , RS (300) ,SFS ( 300 ), THETA $( 30 
10),US( 300 »PS( 300) , RHOS ( 300) ,MS(300) ,VS(300) 

XS(1)=XINITAL 

IF ( IT.GT.I) GO TO 1 

RKSX ( I ) =R KS 

RKS 1 =RKS 

SDERRX ( 1) =SDRX I NT 

N T3=200 

GO TC 2 

RK S i-RKSC 

RKSX ii ) *RKS I 

NT3=NC 

DQ 3 1=2, NT3 
J = I-1 

XS( I )= XS ( J ) +DEL X 

RKSX(I )=RKS1+0.5*SDERRX( IT )*XS(I>**2 

CONTINUE 

RETURN 

END 


C i 
C 2 
C 3 
C 4 
C 3 
C 6 
C 7 
C 8 
C 9 
C 10 
C 11 
C 12 
C 13 
C 14 
C 15 
C 16 
C 17 
C 18 
C 19 
C 20 
C 21 
C 22 
C 23 
C 24 
C 25 
C 26 
C 27 
C 28 
C 29 
C 30- 


SU6R0UTINE SONIC 0 1 

D 2 

THIS ROUTINE COMPUTES SHOCK PROPERTIES IN SUBSONIC REGION UNTIL D 3 

LIMITING SHOCK MACH NUMBER IS COMPUTED. 0 4 

PS=SHOCK PRESSURE, US=SHCCK VELQCI TY( TANGENTIAL )»V S=SHOC K VELOCITY! D 5 

NORMAL ,RHOS= SHOCK DENSITY, MS- SHOCK MACH NUMBER. D 6 

D 7 

REAL MINT , MS IN2 • MS , MSG , MSI D 8 

COMMON /BLK 1/ GAM, MINF, X IN ITAL,2 I NI TAL tRINITAL ,THET AIN, SF IMTL ,L , N 0 9 


It SORXINT.KKUN 

. ABAf ABB, PI f HPA f HP B 

D 

10 

COMMON 

/BLK2/ 

IT f DELX t IBEGt EMtNC t NT 1 f NX tKENDf rtDf PST «UR ( 8) tf 1200) , I 

D 

11 

IP 



D 

12 

COMMON 

/BLK3/ 

XSP( 7) f XT .XTCt SCP fRKS.RKSCtRKSL f SOERRX (ttl ,HSGtXSl 

0 

13 

COMMON 

/BLK*/ 

XSI300) t RKSX (300) »ZS( 300) tRS( 300), SFS( 300) t THETAS! 30 

0 

14 
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10) * US( 300) * P S( 300 ) *RH0S ( 300) » MS (300) *VS(300) D 15 

BL=i(GAM*l. )*MINF>**2 0 16 

DL=4.*GAM*(MlNF**2)-2.*(GAM-l. ) 0 17 

BLD=(BL/CU**(CAR/ (GAM-1. )) 0 18 

THETAMU-THETAIN 0 19 

PRAT=2.0/(GAM*1.0> 0 20 

RAT*(GAM-i.0)/iGAM+1.0) D 21 

GAMM2=1.0/(GAM*MINF**2) 0 22 

IF ( IT.GT.l > GC TO 7 D 23 

UO 1 1 = 1* 7 0 24 

F(I )*1.0/RKSXm 0 25 

1 CONTINUE 0 26 

CALL INT6 (7.0ELX,F.B) 0 27 

NN= 1 0 28 

PM=7 0 29 

GO TO 4 D 30 

2 NN=MM+1 D 31 

MM=MM*2 D 32 

DO 3 1 = 1 * HP 0 33 

P(I)=i.0/KKSX( I) D 34 

3 CONT INUE D 35 

CALL INT6 (MM*DELX,F,B) 0 36 

4 DU 5 I =NN * R M D 37 

THETAS( I >=THETAS( 1 )-F( I) 0 38 

USII )=COS(THETAS(l) ) 0 39 

usm=o.o o 40 

MSIN2=(MINF*SIN(THETAS(I)) >**2 0 41 

PS (I 1=1 PRAT* (SIN (THET AS (I)) ) **2 )- GAM M2 *RAT D 42 

RHOS (I )= ( GAR*1 .0 )/ ((GAM-1.0)*(2.0/MSIN2) ) Q 43 

VSU > = ( 1. /RHOS (I ))*SIN (THETAS U>) D 44 

AS = SORT ( GAM *PS ( 13/ RHO S ( I ) ) 0 45 

MSI I J=(SQRT((US( I)**2)*VS(I)**2))/AS 0 46 

IF (MS(I ) .GT.l.ii GO TO 6 0 47 

5 CONTINUE D 48 

GO TO 2 0 49 

6 NC*I 0 50 

GO TO 10 D 51 

7 DO 8 I *1 * NC D 52 

F(U = l.O/RKSX(n D 53 

8 CONTINUE D 54 

CALL INTfc (NC*DELX*F*8) 0 55 

DO 9 1 = 1* NC 0 56 

THETAS! I)=THETA$(1)-F( I) 0 57 

US i I ) = COS ( ThET AS ( I ) ) D 58 

US(1)=0.0 0 59 

MS I N2= (MINF*SIN( THETAS ( I)))**2 0 60 

PS (I )* (PRAT* (SIN (THET AS ( 1 )) ) **2 ) -GAM R2 *R AT 0 61 

RHOS! I )=( GAMU .0 ) / ( ( G AM- 1. 0 )* ( 2. 0/M S IN 2)) 0 62 

V S ( I ) = ( i . /RHOS ( I ) )*$i N( THE FASH)) D 6 3 

AS* SORT (GAM*PS( I )/RHUS(I ) ) D 64 

MS II ) = (SURT((US( I )**2 )*VS(I >**2))/A$ D 65 

9 CONTINUE 0 66 

PST* PS ( 1 ) *8 LD D 67 

D 68 

PST IS bODT STAGNATION POINT PRESSURE. 0 69 

D 70 

10 RETURN D 71 

END D 72- 
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SUBROUTINE INT6 4NT,OEL,FfB7 

ROUTINE IS 6 PCINT INTEGRATION PROCEDURE. 

DIME NS I UN X 1 200 ) * F4200I 

X(17*0.0 

NN2=NT— 2 

NNl*NT-i 

N=2 

X(27*X( 17*4 (DEL/ 1440. 7*4 (493.*F(N-17 7* 1 133 7. *F I N) 7 - 1 61 8.*F ( N*1 II * ( 
1302.*F(N*2) 7-( 83 .♦FI N*3 ) ) ♦ 19 .*F (N*4 7 7 ) 7 
N=3 

X437*X(27*( (DEL/ I29600.l*( (-4050. *F ( N-2 7 7 * ( t 5430. *F (N- li J ♦ ( 75780.* 
IF (N i 7-4 7G20.*F(N*17 7-4 1170. *F4N*27 7 *4630. *F(N*37 ) 7 1 
DO 1 N*4 »NN2* 1 

X(N7=X(N— 17*4 ( DEL/ 14 400. 7 *( (110.*F4N-377-(930.*F(N-277*48020.*F(N- 
117 7*(8020.*F(N 7 7-4930.*F(N*17 7*41lO.*F(N*27 7 7 7 
CONTINUE 
N=NN1 

X(N7*X(N-1I*( (DEL/ 14400.7*4 (-110. *F( N-4 7 7* ( 770. *F ( N-37 7-4 2580. *F4N 
i— 2 7 7*41 02 20. *F 4 N- 17 ) *4 6370 . *F ( N 7 7 - 4 270. *F 4 N+ 1 7 7 7 I 
N=NT 

X (N7=X(N-17 + 4 ( DEL / 25920. 7*4 4 486.*F4N-5 7 7-4 311 4.*F( N-47 7*4 86 76. *F(N 
l— 377-414.* 64 .*F ( N-2 77*425 68 6.*F4N- 17 7 * 4 85 50 . *F ( N7 7 7 7 
B=X ( NT 7 
OG 2 I *1 » NT * 1 
F l 1 7 =X ( I 7 
CONTINUE 
RETURN 
END 


E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

e 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 
29 
30- 


1 


2 


SUBROUTINE SUBGEN 

THIS ROUTINE COMPUTES SHOCK GEOHE T RY , S TREAM FUNCTION AT SHOCK* AND 
DIVIDES SUBSONIC REGION INTO EQUAL S EGMENT S 4 X SP 4 1 7 7 

RS*RADI AL COMPONENT, 2S*AXIAL COMPONENT *SFS=STKEAM FUNCTION AT SHOC 
K. ALL LENGTHS ARE MEASURED FROM SHOCK. 

COMMON /to LK 1/ G A r* , MI NF , X I NI T AL , 21 N I T AL , R I N I T AL , T HET AI N , SF l N I T L » L , N 
1 ,SDRXINT *RK8N» AB A , ABB * PI * HPA, HP B 

COMMON / BLK2/ I T * DELX , I BEG , EM , NC , NT 1 ,NX * RENC * NO , PST , OR ( 87 ,M2007 , I 
IP 

COMMON / BLK 3/ XSP ( 7 7 » XT * XTC , SCP *RKS * RKSC * R KS1 * SDERRX ( 8 7 * H SC * X S 1 
COMMON / 6LK4/ XS (300 7 * RKSX (300 7 ,ZS( 3 00 7 * RSI 300 7 * SFS( 300 7 • THETAS ( 30 
107 *U SC 300 7 * PS ( 300 7 ,RH0S(30O7 * MS (300 7 *VS(300 7 
RS ( 1 7*R IN IT AL 
2S4 1 7 *21 M T AL 
XS 4 1 7 = X IN I T AL 
SF S ( 1 7 = SF I N l TL 
DO 1 1*1, NC 
F(I7*COS(THETAS(I7 7 
CONTINUE 

CALL I NT 6 (NC*UELX*F,3) 

00 2 1=2, NC 
ZS4I 7=2S( 17*F( 17 
CONT INUE 
DO 3 1 * 1 * NC 
F4 I7*SIN4 THETASU > ) 


F 1 
F 2 
F 3 
F 4 
F 5 
F 6 
F 7 
F B 
F 9 
F 10 
F 11 
F 12 
F 13 
F 14 
F 15 
F 16 
F 17 
F 18 
F 19 
F 20 
F 21 
F 22 
F 23 
F 24 
F 25 
F 26 
F 27 
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3 CONTINUE F 28 

CALL INT6 (NC.DELX.F.B) F 29 

DO 4 1*2, NC F 30 

RS(I)=RS( 1)*F( I) F 3 X 

SFS( I) = ((0.5)**N)*RS( I)**( l*N) F 32 

4 CONTINUE F 33 

IF (IT.&I.l) GO TO 5 F 34 

XT =XS( NC) F 35 

XSl*XT F 36 

DXSP=XT/(NX-1) F 37 

GO TO 6 F 38 

5 XT=XTC F 39 

DXSP*XT/( NX- 1) F 40 

6 DO 7 1=1, NX F 41 

XSPIII =(l-l)*DXSP F 42 

7 CONTINUE F 43 

RETURN F 44 

END F 45- 


SUBROUTINE SUBPRCP G 1 

G 2 

THIS ROUTINE LOOKS UP THE SHOCK PROPERTIES AN U GEOMETRY AT ARC G 3 

LENGTHS EUUAL TO XSP(l) ANU CALLS ROUTINES FRESLOC AND PROPLOC G 4 

G 5 

REAL M INF,MSIN2*MS,MSU,MSI G 6 

COMMON SULK 1/ GAM,MINF,XINITAL,ZINITAL,RIMTAL,THETAIN,SFINITL.L,N G 7 
1 » SDRXINT* RKBN* 7fc A » ABB ,P 1 ,HPA ,HPB G 8 

COMMON /BLK2/ IT , CELX , I BEG , EM , NC, NT l , NX, REND, NO, P ST , OR ( 61 , F ( 200 1 , I G 9 
IP G 10 

COMMON /BLK3/ XS F ( 7 » ,XT ,XT C, SCP , R KS , PK SC , RKSi , SDERRXI d I , HSO , XS 1 G 11 

COMMON /bLK4/ XS 1 300) , RKSXi 300) ,2 S ( 3 00 ) , R S (300 ) , SFS ( 300 I , TH ET AS ( 30 G 12 
10) ,US(300),PS(300),RHOS(300),MS(300I ,VS(300) G 13 

COMMON /BLK6/ SFO ( 7) ,R0 ( 7 ) ,20( 7) , PC ( 7 ) , UO (7 ) , THET AO ( 7 ) , RHOSC ( 7 ) , MS G 14 
10(7),RKQ(7),DELSY (7),RKOC(7), ROCi 7 ) , ZQC 1 7 ) , YO ( 7) G 15 

COMMON /BLK7/ S Y ( 7 , 1 5 ) ,P l ( 7 , 15 ) ,RHi ( 7 , 15 ) , U I ( 1 , 15 ) , MS I ( 7, 15 ) , V 11 7 , G 16 
115) G 17 

DO 1 1=1, NX G 18 

CALL FTLUP ( XSP ( I ) , SFO ( I ) , 1, NC , XS , SF S ) G 19 

CALL FTLLP (XSP1 I ) ,R0 (1) , 1 , NC , XS, RS ) G 20 

CALL FTLUP (XSP(l)»ZU(I),l»NC»XS*ZS) G 21 

CALL FTLUP (XSP(I) ,PQ(I) ,1, NC.XS.PS) G 22 

CALL FTLUP ( XSP ( I ) , VO (I) , 1 , NC , XS, V S ) G 23 

CALL FTLUP (XSP(I),UO(I),l,NC,XS,'jS) G 24 

CALL FTLUP ( XSP( I ) ,THET AO (I) , 1, NC , XS , T HE TAS ) G 25 

CALL FTLUP (XS P ( I ) ,RHOSO (I) , 1 »NC, XS , RHOS ) G 26 

CALL FTLUP (XSP(I) «MSU(I)»l»NC»xS»MS) G 27 

CALL FTLUP (XSF ( I ) , RKO ( I ) , 1 , NC , XS , RK SX ) G 28 

CONTINUE G 29 

CALL PRESLCC G 30 

CALL PRCPLCC G 31 

RETURN G 32 

END G 33- 


SUBROUTINE PRESLCC H 1 

H 2 

THIS ROUTINE COMPUTES THE PRESSUR E( AT EACH XSP VALUE) ALONG RAYS H 3 

NORMAL TO SHOCK. THE PRESSURE EQUATION IS AN APPROXIMATION DEVELOP H 4 

ED BY S.H.MASLEN ACCOUNTING FOR THE CENTRIFUGAL AND RECCMPRESS ICN H 5 

EFFECTS. H .6 

H 7 
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COMMON /BLK 1/ GAM , MI NF , XI NI TAL , ZI N I T AL , R I NI T AL , THET A IN , SF IN ITL , L , N H 8 

1 .SORXINT, RKBN, ABA, ABB, PI ,FPA, HPB H 9 

CGMMON /BLK2/ I T . GEL X , I BEG , EM ,NC . NT 1 ,NX , RENC , NO, PST , OR < 8) , F 1 200 » , I H 10 
IP H 11 

COMMON /BLK5/ RI 1 7 , 1 5 ) , ZI I It 15) . Y I C 7 , 15 ) ,RI C 1 7 , 15 > , Z IC l 7, IS ) , Y IC l 7 H 12 

1,15),RW4C,ZST, E4(8),THB1 50) H 13 

COMMON /BLK6/ SFC 1 7) , RO ( 7) , Z0( 7 1 , PO ( 7 ) , UO (7 1 , 1 HET AO ( 7) » RHCSO( 7 » , M S H 14 
10(7) , RKO ( 7),DELSY(7) »RKOC (7) »ROC( 7),ZUC(7),YQ17) H 15 

COMMON /BLK7/ SY ( 7 ,15 1 , PH 7 , 15 ) , RHII 7 , 15 ) ,U I I 7 , 1 5) , MSI! 7, 15 ) , V I ( 7 . H 16 
115) H 17 

DO 2 1*2. NX H 18 

DELSY( li-SFOm/NTl H 19 

00 1 J = 1 . NT 1 H 20 

SY(I,JI*SFO(I)-J*OELSY(II H 21 

pi(i,j)=pom + (ucm / irkoi n*(RO( i )**n) » > *< sy i i. ji-sfoi u i-usfou h 22 

1)*V0(I )*TAN( THETAOU) )/(RO(I )**N) )*( (1 ./RK0U ))MC0S(THETAQ(I))/R0 H 23 
2( I)) ) )*( ((SYU. J)**2)-SF0(1 1**2)/ (2.*l SF011 1 ** 21 1 1 H 24 

1 CONTINOE H 25 

2 CONT INOE H 26 

RETORN • H 27 

END H 28- 


SOBROOT 1 NE PROFLOC I 1 

I 2 

THIS KOOTINE CCPFOTES TFE PROPEKT I ES ( 0 ENS IT Y= RHI , TANGEN T I AL VELOCI I 3 
TY*UI, NORMAL V EL QC I T Y * V I » AND MACH NUMBER = MSI ) ALONG RAYS NORMAL TO I 4 
SHOCK AT XSP VALUES. AN APPROX IMAT I CN TC NORMAL VELOCITY IS USED. I 5 

I 6 

REAL MINF,MSIN2.MS,MS0,MSI I 7 

COMMON /B LK 1/ G AM , MI N F , X I NIT AL , Z I NI T AC , R I M T AL , THET A IN , SF IN IT L , L , N I 8 
l.SDRXINT.RKBN, AB A , ABB . P 1 ,HP A ,HPB I 9 

COMMON /BLK2/ I T .DEL X , I BEG, EM, NC . NT 1 .NX , RENC , NO, PST , UR ( 8 ) , F ( 200 ) , I I 10 
IP 111 

COMMON /BLK3/ XS P l 7 ) , XT ,XTC , SCP , RKS , RKSC , RKS1 , SDERRX 1 8 ) , HSO ,XS 1 I 12 

COMMON /BLK4/ XS I 300) , RK SX (300) , Z SI 3 00 ) ,RS 1 300) , SFS ( 300 ) , TH ET AS 1 30 I 13 
10),US(300)«PS(300) , RHOS ( 300 ) , MS (300 ) , V S ( 300) I 14 

COMMON / B L K 5/ R I { 7 , 15 ) , Z I ( 7, 1 5) , Y I ( 7 , 15 ) , RI C ( 7 , 1 5 ) , ZI C ( 7 , 15 ) , Y IC ( 7 l 15 

1,15) »RW4C»ZST,E4(a),THB(50) I 16 

COMMON /BLK6/ SFC < 7 ) , RO ( 7 ) , Z0( 7 ) , PO ( 7) ,U0 (7 ) , THET A0( 7 ) , RHCSC ( 7 ) , MS I 17 
101 7) ,RK0(7) , DELSY (7) , RKOC (7),R0C( 7), ZOC( 7),V0( 7) I 18 

COMMON /BLK 7/ S Y ( 7 , 1 5 ) , P I ( 7 , 15 ) ,RHI ( 7 , 15 ) ,U I ( 7 , 1 5 ) , MS I ( 7 , 15 ) , V I ( 7 , I 19 
115) l 20 

PRIT=2.0*GAM/IGAM-1.0) I 21 

HS0=(4.*GAM)/( (GAM + 1. )**2)«-( ( (GAM- 1 • ) / I GAM+1 . ) )**2) *2. / ( (GAM- 1. )*M I 22 
1INF**2) I 23 

DO 3 1=2, NX I 24 

DO 2 J= 1 , NT 1 I 25 

w = SY ( I » J ) l 26 

CALL FTLUP ( to , BP , 2 , NC , S FS , RHCS ) I 27 

CALL FTLUP (W, CP, 2, NC.SFS.PS) I 28 

RHI (I,J)=BP*(PI( I,J)/CP)**(1.0/GAM) I 29 

PRHE=PI(I,J)/RHI(I,J) I 30 

IF (IP. EL. 2) GC TG 1 I 31 

VI( 1 ,J)=VG( I l*(t»/SFC'( i) ) I 32 

UI (1 ,J)=SURT (HSC-FRIT *PRHE-VI (I ,J )**2) I 33 

GO TO 2 I 34 

1 VI(I ,J)=V0<1 )*U1.-(YI (I,J)/YI( I.NT1 )) )**1.1)*( l./i l.-( VI (I ,J)/RKG l 35 

1(1)))) I 36 

UI(I,J)=SCRT(HSC‘-PRIT*PRPE-VI(I,J )**2) I 37 

AE=SORT(GAM*PI ( I , JJ/RHI ( I , J) ) I 38 
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MSI ( I, J) = (SQRT ((bit I ff J J*'*2)+VI< It J)**2) l/AE 

CONTINUE 

CUNT INUE 

RETURN 

END 


I 39 
I 40 
I 41 
I 42 
I 43- 


SUBROUTINE GIP 

THIS ROUTINE COMPUTES THE FLOW FIELD C CORD I NA TES . Y l ^NORMAL CUMPCNE 
NT;RI=RACIAL CCMPCNENT;2I=AXIL COMPONENT. ALL LENGTHS MEASURED FROM 
SHOCK. 

COMMON / 8 LK 1 / GAM , MI NF » XI Ni T AL , ZI N l T AL t R l M T A L » THET A IN , $F I NI T L t L # N 
1 • SDRXINT »RK BN tABAf A 8 BfPIfH PA 9 HP B 

COMMON / BLK 2/ I T » OEL X , I BEG , E M , NC ♦ NT 1 » NX , RENL t ND , P S T , DR ( 8 ) tf (2001* I 
IP 

COMMON /BLK 5/ RI(7fi5)fZI(7,15),YI(7tl5),RlC(7,15),ZIC(7tl5l,YlC<7 
1 1 1 5 ) * Rw4C tZSTf£4(Bi*THB( 50) 

COMMON /B LK 6 / SFC ( 7) » RQ ( 7) *ZO(7)*PC(7J t UC ( 7 ) , ThET AO ( 7 1 * RhGSC( 7 ) * MS 
101 7 ) v RKO ( , DELS Y ( 7) ,RKOC (7) , ROC( 7) ,ZOC( 7) . VO t 7) 

COMMON /BLK 7/ S Y ( 7 , 1 3 ) t P I ( 7, I 5 ) , RH I ( 7 , lb ) t U I ( 7 , 15 J , MS II 7t 15 ) , V I l 1 * 
lib) 

DO 5 I -2 « NX 
S=— DEL SY ( I ) 

F(n*i.o/(RHoscm*uom ) 

DO 1 J- 1 • NT 1 
K = J+ 1 

F(K)*1.0/(RHI( I • J ) *U I ( I v J ) ) 

CONTINUE 

NT2*NTi*l 

CALL INT 6 (NT2*S*F,Bi 
DO 4 J= 1 . NT 1 
K«U J 

IF (N.EQ.O) GO TC 2 

YI( I tJ)*(RO(I)/CCS(THETAO( I) ) ) *(1 . 0- SUR T I i.CM2.0*C0S(THETAC(I))/R 
10(I)**2)*F(K)) ) 

GO TO 3 

YI (I , J)=-F(K) 

Kill * J ) =RQ ( I )- YI ( I * J ) *CCS (THET A0( I ) ) 

RP=RO ( I )-R I ( It J) 

SIGMA* 1.570 79 6326 8- THET AC ( I ) 

Zl =RP/ (TAN( SIGMA)) 

zi(i,j)*zam+zi 

CCNTINUE 

CONTINUE 

RETURN 

ENO 


J 1 
J 2 
J 3 
J 4 
J 5 
J 6 
J 7 
J 8 
J 9 
J 10 
J 11 
J 12 
J 13 
J 14 
J 15 
J 16 
J 17 
J 18 
J 19 
J 20 
J 21 
J 22 
J 23 
J 24 
J 25 
J 26 
J 27 
J 28 
J 29 
J 30 
J 31 
J 32 
J 33 
J 34 
J 35 
J 36 
J 37 
J 38 
J 39 
J 40 
J 4 1- 


SUBROUTINE BAOCURV K 1 

K 2 

THIS SUBROUTINE COMPUTES THE SCALING PARAMETER FOR THE FLCW BY K 3 

K 4 

RATI01NG THE PERTINENT COMPUTED OIMENSICN TC THE KNOWN DIMENSION K 5 

COMMON / B LK 1/ GA M * M I NF , X 1 N I T AL t ZI NI T AL . R I NI T AL , THET A IN , SF I N I TL • L , N K 6 

1* SDRXINT*RKBN, AEA, ABB , PI , HPA ,HPB K 7 
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5 


CGMMCN / 6 LK2/ I T , GELX , I 8EG , EM, NC, NT i ,NX, RENC, NO , PST , OR ( U , F C 200) , I K. 
IP K 

COMMON /8LK3/ XS P( 7) , XT, XTC, SCP,RKS, RKSC, RKS1 , SDERRXt 8 ) ,HSC ,XS 1 K 

COMMON /BLK 5/ RI C 1 ,i 5 ) , Z II 7 , 1 5 I , Y I ( 7 , 1 5 ) ,R I C i 7 , 15 > , ZIC 4 7 , 15 ) , Y IC l 7 K 
1,15) ,KW4C,ZST,E4(8), TH8( 50) K 

K 

IF NX- 4 CCMROL POINTS 2 ANO 3 MUST BE LSED.IF NX*7, CONTROL POINT K 

S 3 AND 5 RECOMMENDED K 

K 

K=3 K 

J* 5 K 

GO TO (3,2,1, 4), L K 

R2=(RI(K,NT1)**2)-RI ( J,NT1)**2 K 

Z2*ZI(K,NT1)-ZI(J,NT1J K 

P*0. 5* (R2/Z 2) K 

K 

PI IS THE KNOWN DISTANCE FROM FOCUS TC DIRECTRIX OF PARABOLA. K 

K 

SCP*P/Pi K 

GO TO 5 K 

C*(RI< J,NT1)**2)-RI (K,NTi)**2 K 

B*(ZI(J,NT I ) **2 )-Z I( K , NT I)**2 K 

A=2.*IZI(J,NT1)-ZI ( K , NT 1 ) ) K 

E*B/ A K 

D-C/ A K 

G-E-ZI ( K* NT 1 ) K 

Cl*( AB8**4)/D**2 K 

C2-{ 2./A6B**2) *G*D+I ( RI ( K ,NT 1) **2 ) / ( AB8**2) 1-1. K 

C3*C 1*C2 K 

C4=( ABB**4)*(<G/C)**2) K 

AX* ( -C3+SUR Tl(C3**2)-4.*C4) 1/2. K 

ABC* SQRT ( AX ) K 

K 

ABA IS AXIS OF KhCwN ELLIPSE PARALLEL TC FLCw K 

K 


8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 


SCP* ABC/A BA 
GO TO 5 

RB2»(RUK,NT1)**2J*ZI U,NT1)**2 
RB3=(RI( J,NT1 )**2)*ZI ( J , NT i ) **Z 
ZB2*Z I (K,NT ll-ZI ( J , N T1 ) 

RB23*R82-KB3 
ZN*RB23/(2.0*ZB2 J 

KKBC=SmRT( ( ( Z 1 (K,NTli-ZN )*** J+RI 4 K »NT 1 1**2) 

RKBN IS THE KNCUN RADIOS OF CURVATURE CF SPHERE 

SCP*RKBC/RKBN 
GO TC 5 

B*(RI(K,NT1 ) **2)-R 1 1 J ,NT 1) *%2 
C=(ZIl J,NT1 )**2)-ZI(K,NTl)**2 
0=2. *( ZI(K,NT1)-ZI CJ,NTi) ) 

E=b/D 

H*C/D 

G*ZI (K , NT 1 ) ♦H 
Cl=( HP B**4)/E**2 

C2* (2./HPB**2 J *G*E-( (RI(K,NT i)**2)/HPB**2)-l. 
C3*C 1*C2 

C4*(HPB**4)*UG/E)**2) 

HX=(-C3*S4RT((C3**2)-4.*C4n/2. 

HY*SURT(HX) 

HP A IS AXIS OF KNC*N HYPERBOLA PARALLEL TC FLOW 

SCP* HY/HP A 

RETURN 

END 


K 43 
K 44 
K 45 
K 46 
K 47 
K 48 
K 49 
K 50 
K hi 
K 52 
K 53 
K 54 
K 55 
K 56 
K 57 
K 56 
K 59 
K 60 
K 61 
K 62 
K 63 
K 64 
K 65 
K 66 
K 67 
K 68 
K 69 
K 70 
K 71 
K 72 
K 7 3- 
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SU6R0UT INE SCALE L 1 

L 2 

THIS ROUTINE MODIFIES ALL LINEAR DIMENSIONS wITH SCALING PARAMETER L 3 
(SCP). L 4 

L 5 

COMMON / BLK 1/ GAM,MINF,XINITAL»ZINITAL,RINITAL*THETAIN»SFINITL*L»N L 6 
1,SDRXINT,RKBN,AEA,ABB,P1,HPA,HPB L 7 

COMMON /BLK2/ IT , 0 ELX , I B EG , EM , NC. NT 1 , N X.RENC , NO ,PST ,DR ( 6) , F I 200) . I L B 

IP L 9 

COMMON /6LK3/ XS P I 7 I , XT , XTC, SCP * RKS . RKSC , RKS 1 , SDERRXi B I ,H SO « XS l L 10 

COMMON /BLK4/ X$ l 3 00 ) , RKSX (300) ,ZS ( 300 ) ,RS 4 300 ) , SFS ( 300 ) , TH ET AS ( 30 L 11 
10),US(300),PS(300) »RHOS (300),MS(300)»VS(300) L 12 

COMMON /BLK5/ R I ( 7 , 1 5 1 , 21 ( 7 , 1 5 1 , Y I < 7 , 1 5 ) , R I C I 7 , IS i , L IC l 7 , 1 5 ) , TIC ( 7 L 13 

1>15) ,Rr4C »ZST»E4(B)*TH8(50) L 14 

COMMON /BLK6/ SF C I 7) , RO < 7 ) , ZC( 7 1 , PC I 71 . UO < 7 I , THET AO ( 7 ) , RHCSOI 7 1 , MS L 15 
10(71, RKO(7),OELSY( 7) , RKOC ( 7 ) ,RUC( 7) ,ZGC< 7) ,V0(7) L 16 

IF (L.Ew.2.0R.L.EC.4> GO TO 1 L 17 

SCP a 1.0/SCP L 18 

RKSC a RKSl* SC P L 19 

DO 2 l a 2»NX L 20 

RIC( I »NT1)*RI( I, NT 1 ) *SC P L 21 

ZIC(I,NTl) a ZI( I»NT1I *SCP L 22 

YIC(I»NTl) a YI(I»NTl)*SCP L 23 

zocm*zcm*scp l 24 

ROC ( l)=RO(l)*SCP L 25 

RKOC(I) a RKC(I)*SCP L 26 

CONTINUE L 27 

DO 3 1 a 1 » NC L 28 

xsi n*xs( ii*scp l 2 9 

RSU ) a RS(I)*SCP L 30 

2S(I ) a ZS( I )*SCP L 31 

RKSX (I ) a RKSX(I )*SCP L 32 

CONTINUE L 33 

RETURN C 34 

END L 35- 


SUBROUTINE SOSCD M 1 

M 2 

THIS ROUTINE COMPUTES THE STAGNAT ICN-PCINT SHOCK STANO-OFF M 3 

DISTANCE (ZST)ANC THE RESPECTIVE BODY ANGLES! THB ) . THB IS ANGLE OF M 4 

INCLINATION OF BOOY SURFACE RITH LONGITUDINAL AXIS. M 5 

M 6 

COMMON /B LK 1/ GAM , MI NF , X IN I T AL, i I N I T AL , R IN I T AL ,T HET AIN, SF I N ITL , L , N M 7 

1,SDRXINT,RK6N, AB A , ABB , P 1 , HPA , HPB M 8 

COMMON /BLK2/ I T ,DELX , I BEG , EM, NC, NT 1 ,NX , RENC, NO, PST , UR ( 8) ♦ F ( 2001 , I M 9 

IP M 10 

COMMON /BLK5/ R1 ( 7 , 15 ) , 21 ( 7 , 1 5 1 , Y 1 1 7 , 1 5 i , RI Ct 7 , 15 J , L IC ( 7 , 15 I , Y IC ( 7 M 11 

1,15),RR4C,ZST,E4(6),THB(50) M 12 

M 13 

IF NX a 4, CONTROL POINT 2 MUST BE USEO.IF NX a 7,CCNTR0L POINT 3 RECC- M 14 

MENDED. M 15 

M 16 

K=3 M 17 

GO TO 13,2,1,4), L M 18 

1 Bl a ZIC(K,NTl) M 19 

B2=-((RIC(K,NT1)**2)/(2.*P1)) M 20 

GO TO 5 M 21 

2 Bl a ZIC(K*NTl)— AEA M 22 

B2=SQRTUABA**2)-UABA/ABB)**2)*RIC(K,NT1)**2) M 23 

GO TO 5 M 24 

3 81 a Z IC ( K» NT 1)-RK8N M 25 
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B2*$WRT( IRKBN**2)-RIC1K,NT1)**2) 

GO TO 5 

4 B1*ZIC ( K, NT 1 ) +HP A 

82 s — SURT (1HPA**2)+((HPA/HPB)**2I*RIC 1K,NT1)**2) 

5 ZST s BU-B2 
LM=NX-1 

DO 10 I s ! f LH 
J=Ui 

GO TO 16.7,8,9), L 

6 THB1I) *ATAM 1 RKBN- 1 2 1 C ( J , NTi ) -ZST ) ) / R I C ( J , NT 1 ) ) 

GO TO 10 

T THB 1 1 ) =AT AM ( 1 ABB/ A8A)**2) * ( A8A-Z 10 ( J, Ml) + 2ST)/RIC(J»NUiJ 

GO TO 10 

8 THB(I) S ATAN(P1/RIC(J,NT1)) 

GO TO 10 

9 THB 1 1 ) *AT AN 1 1 1 EPB/HP A )**2 ) *1 HPA+Z 1C 1 J, NT 1WST )/R ICC JvNTl )) 

10 CONTINUE 
RETURN 
END 


M 26 
M 27 
M 28 
M 29 
M 30 
M 31 
M 32 
M 33 
M 34 
M 35 
M 36 
M 37 
M 38 
M 39 
M 40 
M 41 
M 42 
M 43 
M 44- 


SUBROUTINE R4E4 N 

N 

THIS ROUTINE COMPUTES THE ANALYTICAL 8CCY RADIUS AT LAST BODY PO IN N 
T ( NX ,NT 1 ) WITH CCMPUTED ZIC(NX,NTi> ANC ZST. N 

ALSO ERROR BETWEEN CCMPUTED BODY 1R IC 1 NX ,NT 1) ) AND ANALYTICAL N 

BODY 1 RW4C ) I S DETERMINED. N 

N 

COMMON /bLK 1/ GAM, MINE, XINI TAL , Z l NI T AL ,RI M TAL , THE TA I N, SE I NITL ,L ,N N 
1 , SDRXINT , RKBN, ABA , ABB, PI ,hPA , HPB N 

COMMON /BLK2/ 1 T ,DELX, I BEG , EM ,NC , NT 1 , NX, REND, NO , PST ,DR ( 8) , F (200) , I N 
IP N 

COMMON /BLK 5/ RI (7, 15), 2117, 15) »YI( 7 ,15) ,RIC( / ,15) , Z IC17,15),YIC<7 N 
1,15) ,RW4C, ZST,£4 (8),THB1 50) N 

GO TO 13,2, 1,4), L N 

RW4C s SQRT(2.*Pl*(ZIC(NX,NTi)— ZST) ) N 

GO TO 5 N 

63 s S wRT 1 1 ABA* *2 ) - < 1 AB A+LST )— ZIC1NX,NT1 i 1**2 J N 

RW4C S ( ABb/ABA)*83 N 

GO TO 5 N 

b3 s (ZlCtNX,NTli-(RK8N+Z$T))**2 N 

RW4C=S WRT (<RK6N**2)-B3) N 

GO TO 5 N 

B3 s SgRT 1 ( UICtNX,NTll*HPA-ZST)**2)-HPA**2) N 

RW4C- ( HPfi /h PA ) *£ 3 N 

E41 IT)*t l R I C (NX , NT 1 ) /RW4C ) — 1 • ) *100. N 

RETURN N 

END N 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 
27- 


SUBROUTINE CORRI 0 1 

0 2 

THIS ROUTINE IS USEC TO CORRECT THE LENGTH ALONG SHOCK TO SONIC 0 3 

REG I ON IX TC) If THE VALUE CF SDERRX IS POSITIVE FOR THE PARTICULAR C 4 

I TER AT ION { IT ) . 0 5 

0 6 

COMMON / BLK2/ IT , DELX , I BEG, EM , NC , NT 1 ,N X , RENC , NO , P ST , DR ( t ) , F ( 2C0 ) , t 0 7 
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ip □ a 

COMMON / BLK3/ XS P ( 7) , XT , XTC, SCP,RKS, RK SC f RKS 1 , SDERRX( 8 ) »H SO , XS1 C 9 

SDERRX2=0.5*S0ERPX(I T ) 0 10 

R I G= SUERRX2 /RK SC 0 11 

SRI G = SwRT (RIG) 0 12 

XRK= (XS1*RKSC)/RKS 0 13 

XTC=SQRT( 1. 0/RIG)* TAN ( SRI G*XRK) 0 14 

DELX=XTC/(NC-1) 0 15 

RETURN 0 16 

END C 17- 


SUBROUTINE CQRR2 P 1 

P 2 

THIS ROUTINE IS USED TO CORRECT THE LENGTH TO THE SONIC POINT(XTC) P 3 

IF SDERRX IS NEGATIVE FOR THE ITERAT IGN( IT) . P 4 

P 5 

COMMON / B LK2/ I T , CELX, I BEGt EM, NC, NT 1 *NX* RENC, ND,PST,DR< £ ) ,F ( 2C0 ) , I P 6 

IP P 7 

COMMON /BLK3/ XSP ( 7 ) , XT ,XTC» SCP,RKS t RK SC ,RKS1 , SDERRX! 8) ,H SO ,XSi P 8 

SDERRX2*0.5*SDER RX 4 1 T i P 9 

srkss= sort (absiscerrx2)/rksc) p 10 

SRKSD* SORT (RKSC/ABSISDERRX21 ) P 11 

XRK* (XS1*RKSC)/RKS P 12 

EXP1=EXP(2.0*SRKSS*XRK) P 13 

EXPN*EXP1-1.0 P 14 

EXP0=EXP1 *1 *0 P 15 

EXPU*EXPN/EXPO P 16 

XTC*$RKSC*EXPU P 17 

DELX*XTC/(NC-1) P 18 

RETURN P 19 

END P 20- 


SUBROUTINE SOCORR 0 1 

0 2 

THIS ROUTINE COMPUTES SDERRX IF IT IS GREATER THAN 2 C 3 

0 A 

COMMON /B LK2/ I T f OELX , I BEG , EM , NC , NT 1 #NX »RENC f NO t PST , DR < 8) t F (200 ) * I 0 5 

IP C 6 

COMMON /BLK3/ X$P( 7) , XT , XTC, SCP ,RKS ,RKSC,RK$1 , SOERRXI 8 ) »HSG ,XS1 0 7 

COMMON / 6LK5/ R II 7 # 1 5 1 , 1 I ( 7, 15 ) f V I ( 7 , 1 5 ) ,R I Cl 7 # 15 ) , L IC I 7* 1 5 ) , 1 1 C i 7 C 8 
1 ,15) ,R*4C,ZST,E4(8),THe(50) 0 9 

I T * I T + 1 C 10 

IF (IT. EG. 9) SDERRXI U*SCERRX(7) 0 11 

IF (IT.EC.9) SCERRXI 2 ) a SDERRX I 8 ) C 12 

IF l IT. EG. 9) E 4 ( I ) -E4 (7) 0 13 

IF (IT. EG. 9) E4( 2)=E4( 8) Q 14 

IF (IT.EC.S) I T* 3 Q 15 

SDERRX ( IT ) = SOERRXi IT-1)-E4( IT-1)* ( SO£RRX( IT-1 I-SDERRX ( I T-2 ) ) / ( E4 ( I 0 16 

iT-l)-E4( IT-2)) 0 17 

RETURN 0 18 

END G 19- 
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APPENDIX B - Continued 


SUBROUTINE OUTPUT R 1 

REAL MINE, P S I N2 , MS ,MSO , MSI R 2 

COMMCN /BLK2/ IT , OELX , I BEG, EM.NC , NT 1 , NX, REND, NO, PST , OR I B) ,F ( 200) , I R 3 
IP R A 

COMMON /BCK 3/ XSP ( 7 1 , XT , XTC , SCP ,K KS , RK SC ,RKS1 , SDERKX < 8) ,H SC ,XS1 R 5 

COMMON /6 LK 5/ RI(7 , 15 ) , 21 1 7, 15) ,Y I ( 7 , 15 ) , RI C< 7 , 15 ) , 2 IC i 7, 15 1 , V IC < 7 R 6 
1 » 151 ,Rw4C ,2 ST* E4 ( b ) » THB( 50 ) R 7 

COMMON /BLK6/ SFC ( 7 ) , RG 4 7 ) , 20( 7 ) , PO < 7) , UO (7 ) , T HET AO ( 7 ) , RHOSO( 7 ) , MS R B 
10( 7) ,RKO( 7) ,DELSY( 7) ,R*0C(7) ,ROC( 7) ,20C( 7) , VQ ( 7) R 9 

COMMON /BLK 7/ S Y ( 7 » 1 5 ) , P I (7 , 15 ) ,RHl C 7 , 15 ) ,U I ( 7 , 15 ) , MS I < 7, 15 ) , V l ( 7 , R 10 
115) R 11 

WRITE <6, 7) RKSC.2ST R 12 

WRITE <6,6) PST R 13 

WRITE (6,9) R 14 

DO 1 1=2, NX R 15 

WRITE (6,10) ASP ( 1 ), 2UC ( I ) , RUC ( I ) « 2 IC( I,NT1),RIC(I,NT1) » Y IC ( I *NT 1 ) R 16 
1 , RKOC ( I ) R 17 

1 CONTINUE R 18 

WRITE (6,11) R 19 

DO 2 1=2, NX R 20 

J=I-1 R 21 

WRITE (6,10) MSI(I,NT1),P0U),UCI< I ) , V0( I ) ,RHUSQ< I ) ,MSO(I) ,TF8i J) R 22 

2 CONTINUE R 23 

DO 6 1=2, NX R 24 

WRITE (6,15) XSF(I) R 25 

IF ( I.GT.2) GO TO 3 R 26 

WRITE (6,12) R 27 

3 DO 4 J= 1 , NT 1 R 28 

WRITE (6, 16 ) PI ( I, J), RHI ( I, J ),UI ( I, J), Vi ( I, J) ,MSI II ,J) ,57(1 ,J) R 29 

4 CONTINUE R 30 

WRITE (6,13) R 31 

DO 5 J = 1 , NT 1 R 32 

WRITE (6,14) TI< I, J) ,2I< I, J) ,RI (I ,J) R 33 

CONTINUE R 34 

CONTINUE R 35 

RETURN R 36 

R 37 

FORMAT </3X,22H STAG PT SHK RAD CUR = , El 1 .4 , 5X , 28H STAG PT STAND 0 R 38 
IFF 01ST =,E11.4) R 39 

FORMAT (3X,19HBCDV STAG PT PRESS= , E 1 1 . 4 ) R 40 

FORMAT ( //4X, 10FARC LENGTH ,3X , 1 5HSHCCK AXIAL L OC , 3X , 16HSHCCK RAD IA R 41 
1L L0C,3X, 14HB0CY AXIAL LOC , 3X , 15H80DY RADIAL COC , 3X , 15HB00 Y NORMAL R 42 
2 LOC, 3X, 2 Oh SHOCK RADIUS OF CURV) A 43 

10 FORMAT (3XElU4,5XEll.4,7XE11.4,7XE11.4,6XE11.4,7XE11.4,8XEll-4) R 44 

11 FORMAT ( //4X, 10FMACH NO WL ,3X , 1 5HSH0CK PRESSURE , 3X, 16HTAN SHK VELG R 45 

1CITY,3X,16FN0R SFK VE LOC I T Y , 3X , 14HSHCCK DENSITY, 3X, 15HSF0CK MACH N R 46 
2UMB, 6X, lOHbCOY ANGLE) R 47 

12 FORMAT (//3X,14FF IELD PRESSURE , 4X , 13HF I ELD DENS 1 T Y , 4X, 2CHF I ELD VEL R 48 

1UC I T Y< TANG) ,4X,20HFIELD VE LOC I T Y ( NORM ) ,4X , 13H F I ELD MACH NC,4X,15HF R 49 
2 IELD STRP FUNC) R 50 

13 FORMAT ( //3X,55HFIELU NORMAL LOC FIELO AXIAL LOC FIELC RADI A R 51 

1L LOC) R 52 

14 FORMAT (5XE11.4,6X£11.4 f 6XE11.4) R 53 

15 FORMAT (//3X.5F XSP=,E11.4) R 54 

lo FORMAT (5XE11.4,6XE11.4,7X£11.4,14XE11.4,10XE11.4,7XE11.4) R 55 

END R 56- 
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APPENDIX B - Continued 


SUBRCUT INE SUPER S 1 

S 2 

THIS ROUTINE COMPUTES THE SUPERSONIC PLOW FIELD. S 3 

S 4 

REAL M INF,MSIN2,MS,MS0,MSI S 5 

COMMON /BLK 1/ GAM, M INF, X INIT AC , Zl N IT AL ,R IN I TAL , THETA IN, SF I N I TL ,L , N S 6 

1 , SOR XI NT ,RKBN,ABA,ABB,P1,HPA,HPB S 7 

COMMON /6LK2/ I T , DELX , I BEG, EM,NC, NT I ,NX,RENC, NO, PST ,OR l 8) , F I 2C0) , I S 8 

IP S 9 

COMMON /BLK37 XS M 7 ) , XT , XTC, SCP, RKS , RK SC ,RKS1 , SUERRXl 8 ) ,H SO , XS 1 S 10 

COMMON /BLK4/ XS (300) ,RKSX{ 300) ,2 S (3 00 ) , RS (300) , SFS 1300 ) , TF ET AS 1 30 S 11 

10),US(300),PS(300),RHOS(300),MS(300),VS<300) S 12 

COMMON /bLK 5/ RI ( 7 ,1 5 ) , Z 1 1 7, 1 5 ) ,Y I 1 7 , 1 5 I ,R I C I 7 , 15 ) , Z IC I 7, 15 I , Y I C I 7 S 13 

1,15),RW4C,ZST,E4(8)*THB(50) S 14 

COMMON /BLK6/ SF C I 71 , PO 1 7 » , 20 1 7 I , PC I 7 ) ,GC 1 7 I , THET AO 1 7 1 , RHO SC I 7 ) , MS S 15 

10(7) *RKO( 7),GELSY( 7), RKOC (71, ROC I 7),ZOC(7),VO(7) S 16 

COMMON /BLK7/ ST 1 7 , 1 5) , P 1 1 7 , 15 ) , RHI I 7 , 15 ) ,U 1 1 7 , 15 ) , MS 1 1 7, 1 5 ) , V 1 1 7 , S 17 

115) S 18 

01 MENS ION RBI20), ZBI20) S 19 

I0PT=0 S 20 

THC=0. S 21 

THC=0.0174532S252*THC S 22 

S 23 

IF IOPT=0, ANALYTIC CALCULATION CONTINUES THROUGH SUPERSONIC REGICN S 24 

IF I OPT = 1 » A CONE JUNCTURE WILL BE MATCFEO TC ANALYTIC SHAPE. THC = S 25 

BLUNT CONE HALF ANGLE. IF THC=0 ANO IQPT*1,A CYLINDER WILL BE MATCH S 2t> 

EO TO ANALYTIC SHAPE. S 27 

S 28 

DETHB* .0 S 29 

I TE S T* 1 S 30 

NH=0 S 31 

NHC=0 S 32 

NHH=0 S 33 

THS1=0. S 34 

1 SUM =NX- 1 S 35 

00 1 J*l» ISUM S 36 

1=3*1 S 37 

R6(J)=RIC(I,NT1) S 38 

ZB( J)=*ZIC(I,NTl) S 39 

CONTINUE S 40 

ICQUNT=Q S 41 

RHS=0. S 42 

R 1=0.0 S 43 

THSC=ASIN(SQRT((1./MINF**2)*((GAM*1.)/2.)*SIN(THC)**2)) S 44 

S 45 

THIS EOUATICN COMPUTES SHARP CONE SHOCK ANGLE. S 46 

S 47 

NP=0 S 48 

PRIT=2.0*GAM/(GAM-1.0) S 49 

GAMM2=1.0/(GAM*MINF**2) S 50 

RAT=(GAM-1.0)/(GAM*1.0) S 51 

PRAT *2 .0/ (G AM-t- 1. 0) S 52 

OR ( 1 )=• 1 S 53 

IT *1 S 54 

0ETH=(THETAS(NC)-TFETAS(NC-1) )/1.5 S 55 

I/«0 S 56 

ND=NC*7 S 57 

10=0 S 58 

RKSX ( ND) = RKSX (NC)*CR( IT) S 59 

S 60 

THE SHOCK RAOIUS OF CURVATURE IS COMPUTED S 61 

S 62 

RKC2=0.5*(RKSXIN0)*RKSX(NC) ) S 63 

NA=NC+1 s 

IF (IT.GT.l) THETAS(NO)=THETASINC )*IRKXS/RKC2 ) S 65 

IF (IT.GT.l) OETF=(THETASINO)— THETAS INC) )/7. S 66 
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APPENDIX B - Continued 


DO 4 1*1, 8 s 67 

F1I)*-RKC2 s 68 

4 CONTINUE S 69 

CALL INT6 <8,0ETH,F,B) S 70 

S 71 

THE FOLLOWING CC LOOP COMPUTES SHOCK PPOPEKTIES S 72 

S 73 

DO 5 I=NA,ND S 74 

THETAS! I)=THETAS1I-1)+CETH S 75 

M= I— NC ♦ 1 S 76 

XS1I )*XS1NC)«-F4M) S 77 

IF ( IZ.EU.l.ANC.RKSXlNC) .GT.O.) A S ( I 1 = XS I NC ) -F 1 M) S 78 

OSm*COS!THETAS(I>) S 79 

MSIN2=1MINF*SIMTHETAS1I1) 1**2 S 80 

PS 1 1 I * < PRAT* SI N 1 THET AS (I i 1 4*2 I -GA MM2* PAT S 81 

RHOSII )*4GAM+1. 01/11 GAM-1.01+12.Q/MSIN2) I S 82 

VS(I )*ll.0/RH0S4 I) )*S INITHET AS ( 1 1 ) S 83 

AS* SORT (GAM*PS I I I /RHOS ( I ) ) S 84 

MSIII=4SttRTHuSlI)**2)*VS4I>**2))/AS S 85 

CONTINUE S 86 

IF 4IT.E0.il RKXS=RKC2*4THETAS4ND)-TFiETAS(NC)l S 87 

S 88 

THE FOLLOWING ROUTINE CCKPUTES THE SHOCK GEOMETRY. S 89 

S 90 

DELX=tXS4NC)-XSlNC))/7. S 91 

DO 6 I =NC »ND S 92 

M=I-NC*1 S 93 

F I M I *COS( THETAS! I) I S 94 

CONTINUE S 95 

CALL IN T6 ( 8, DEL X ,F «8 i S 96 

DO 7 I=NA» NO S 97 

M=I-NC*1 S 98 

2S 4 1 ) = ZS INC 1 *F (M 1 S 99 

CONTINUE S 100 

DO 8 I*NC,ND S 101 

M= I-NC + 1 S 102 

F !M)*SIN!TFETAS! II I S 103 

CONTINUE S 104 

CALL INT6 18,DELX,F,B) S 105 

DO 9 I =NA , NO S 106 

M* I- NC ♦ 1 S 107 

RSll)=RS4NC)*F4M) S 108 

SFS4I1*110.5**M*RS4 11**1 1+N) I S 109 

CONTINUE S 110 

S 111 

THE FOLLOWING ROUTINE COMPUTES THE PROPERTIES ALONG A RAY S 112 

S 113 

S = SF S l ND ) / NT 1 S 114 

I TEST* 1 S 115 

10 DU 12 J*1,NT1 S 116 

IF ! I TEST. EC. 2) GO TO 11 S 117 

SY14,J)=SFS1N0)-J*S S 118 

PI 14, Ji =PS!ND) ♦tuSIND) /< RKSX 1N0I* IRS 1NC)**N) 1 I * 1 S Y ( 4, J I -SF S 1NDI I - 4 S 119 

1 ISPS 1ND)*VS 1NU)*TAN1 THET AS 4ND 1 1 / 1 R S 1 NO I **M 1* 4 1 1. /RKSX 4 NO I I ♦CCSIIH S 120 
2ETAS(ND)i/RS4NDl 11*1 1 (SY 14»vil**2)-lSFS 4 NO 1**2) 1 / 12. *1 SFS4 ND 1 **2) 1 1 S 121 
W = S Y 1 4 » J 1 S 122 

CALL FTLUP 1 w , BP , 1 ,NC , SFS , RHOS I S 123 

CALL FTLUP t W, CP , 1 ,NC , SFS ,PS I S 124 

RHI 14, J)=BP*1PI14, J)/CPl**ll. O/GAM) S 125 

PRHE=PI(4,J1/RHI44»J1 S 126 

VI 44,J)*VSlND)*lw/SFS 4ND) I S 127 

UI14,J)*SjRT (HSC-PRI T*PRHE-VI (4,J)**2) S 128 

GO TO 12 S 129 

11 W*S Y 1 4 » J 1 S 130 

VI 44,J)*VS1 N01 *1 I1.-4YI 44, J) /YI44,NT1) 1 1 1* 11./ 11.-4 Y 144, J »/ RKSX 4ND S 131 

111)1 S 132 
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PRhE*PIl4,J)/RHll^.J) S 133 

UI(4r J)*S uRT1hS0-PRIT*PRHE-VI(4,J 1**2) S 134 

AE=SgRT(GAF*PRFE 1 S 135 

MSl<4,Jl*(SURT(<0ll4,J)4*2)*VI(4, J)**2)»/AE S 136 

12 CONTINUE S 137 

C S 138 

C THE FOLLOWING STATEMENTS COMPUTE THE FLO* FIELD GEOMETRY. S 139 

C S 140 

F(1)=1.0/1RP0S(NC1*US tND I ) S 141 

DO 13 J- 1 1 K T 1 S 142 

K=l«-J S 143 

FIK)*1.0/(RHU4,J»*UI (4, Jl) S 144 

13 CONTINUE S 145 

NT 2*NT 1* 1 S 146 

CALL INT6 (NT2,S,F,6» S 147 

00 16 J=1,NT1 S 148 

K= 1 ♦ J S 149 

IF (N.Eu.O) GO TC 14 S 150 

YI (4,J1=IRSIN0)/CGS(THETAS<N01 II* 1 1.0- SORT! 1. 0- 1 2 . 0*C0S I THE T AS (NCI S 151 

1)/RS(N0)**2)*F(KJJ ) S 152 

GO TO 15 S 153 

14 Y I (4 * Jl *F( K) S 154 

15 RI(4,J1=RSIND1-YI(4,J|ACCS(THETAS(NC11 S 155 

SIGMA® THETA IN-THE T AS INCI S 156 

RP=KS (NO) -R I (4, Jl S 157 

Z1*RP/ (TANIS IGMA11 S 158 

ZI(4,J)*ZS(N0)*21 S 159 

16 CONT INUE S 160 

IF I I TEST. EG. 2 ) GO TO 17 S 161 

Y S* Y I ( 4 • NT 1 ) S 162 

I TES T*2 S 163 

GO TO 10 S 164 

C S 165 

C THE FOLLOWING STATEMENTS SCALE THE FLOw FIELD S 166 

C S 167 

17 IF (A8SIYS-YII4,NT1)1.LE.. 000011 GO TO 18 S 168 

Y S*V I ( 4 , NT 1 1 S 169 

GO TO 10 S 170 

18 RIC(4»NT1)*RI(4,NT 11 S 171 

ZIC(4,NT1 1*21(4, NT1) S 172 

YIC(4»NT11*YI(4»NT11 S 173 

ZGC( 41 *ZS INC) S 174 

ROC ( 41 *RS ( NG 1 S 175 

RKOC (4 1*RKC 2 S 176 

C S 177 

C THE FOLLOWING STATEMENTS COMPARE THE CCMPUTEO AND ANALYTICAL S 178 

C BODY LCCATICNS ANC CETERMINE NEW CORRECTION S 179 

C S 180 

IF (NHC.Et.lJ GC TC 34 S 181 

IF (Rl.EU.0.01 GC TC 19 S 182 

IF (R1.GT.0..AND.NP.EC.01 GO TO 30 S 183 

IF INP.EU.l) GC TO 35 S 164 

19 GO TO 120,21,22,231, L S 165 

20 B3*(IZST*RK8M-ZIC(4,NT1)1**2 S 186 

RW4C*SGRT (1RKBN**2)-B3) S 187 

GO TO 24 S 188 

21 63*1 IZST*ABA1-ZICI4,NT1) 1**2 S 189 

RW4C*(A6tt/ABAI*SGRT(( ABA**21-B31 S 190 

GO TO 24 S 191 

22 Rw4C=SgkT(2.*Fl*IZIC(4,NTl)-ZST)) S 152 

GO TO 24 S 193 

23 B3*(ZIC(4,NT1)*(HPA-ZST1 1**2 S 194 

Rw4C*(HP6/FPA)*SCRT(B3-(HPA**211 S 195 

24 £4(IT1=((KIC(4»NT1)/RW4C1-1.)*100. S 196 

IF IABS(E4( mi.GT.EM) GC TO 3o S 197 

1 SUM* I SUM + 1 S 198 

I * I S UM S 199 
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GO TO 427,26,25,281, L S 200 

25 THB( I) 3 ATAN(P1/RIC(4,NT1)) S 201 

GO TC 29 S 202 

26 THB(l) 3 ATAM(( AB8/ABA)**2)*(ABA-ZIC(4,NT1)*ZST)/RIC(4,NT1)) S 203 

GO TO 29 S 204 

27 THB<I)=ATAN((KKBN-(ZIC(4,NT1)-ZST))/RIC(4,NT1)) S205 

GO TO 29 S 206 

26 THBlI) 3 ATAM({HPB/PPA)**2)*(HPA«-ZICI4,NTl)-2ST)/iUC(4,NTl)) S 207 

29 IF (THC.EC.O..ANC.IOPT.EU.I) GO TO 31 S 208 

IF (IOPT.EC.O) GC TO 42 S 209 

RB ( I SUM )=R 1C(4,NT1) S210 

ZBI I$UM)=ZtC(4,NTi) S 211 

IF ITH8U5GM1.GT.THC) GO TO 42 S 212 

CALL FTLUP ( THC , R 1 ,1 , I SUM ,THB, R8) S 213 

CALL FTLUP 1 THC , 22 , 1 , I SUM , THB ,ZB) S 214 

OETC-UETH S 215 

C S 216 

C THE BOOY-CCNE JUNCTURE IS LOCATED. R1 ANC Z2 ARE THE RAOIAL AND AXI S 217 

C AL LOCATIONS R ESP ECT I VELL Y OF THE JUNCTURE POINT. S 218 

C S 219 

WRITE (6,55) S 220 

OETH=OETH*(Rl-RE( ISLM-1) )/ IRBi ISUM)— RBI ISUM— l ) ) S 221 

DR(i)=.5*CR( IT) S 222 

IT«1 S 223 

GO TO 3 S 224 

30 RW4C*R1 S 225 

E4(IT)=((RIC(4,NT1)/RW4C)-1.0)*100. S 226 

IF (ABS(E41IT)).CT.EMi GC TO 36 S 22 7 

NP 3 1 S 228 

l)ETH*UETC S 229 

C S 230 

C THE BODY- CONE JUNCTURE AND RESPECTIVE SHOCK POINT ARE LOCATED. S 231 

C THE PROPERTIES ALONG RAY INTERSECTING JUNCTURE POINT ARE CGMPGTEC. S 232 

C S 233 

GO TO 42 S 234 

31 IF (NH.EU.l) GO TO 32 S 235 

C S 236 

C BOOY-C YL INCER JUNCTURE IS COMPUTE 0. NOT E IN THE FOLLOWING STATE S 237 

C MENTS THE F IRST , S£ VENTH, AND EIGHT SHCULC BE CHANGED IF BODY OTHER S 238 

C THAN A SPHERE IS CONSIDERED. S 239 

C S 240 

IF ( (RK6N*Z$T-ZIC(4,NTl) ).GT.O. ) GO TO 42 S 241 

DETHB»OETH S 242 

NH 3 1 S 243 

OETH 3 . 5*DET H S 244 

DR ( 1 ) 3 . 5*CR ( IT I S 245 

IT*1 S 246 

GO TO 3 S 247 

52 IF ( (RKBN*2ST-ZIC(4,NT1) ).GT..000l) GO TO 42 S 248 

IF (ABS(RK8N<-ZST-ZIC(4,NT1)).LE..0001) GO TC 34 S 249 

IF (NHH.EU.l) GO TO 33 S 250 

DR ( l ) 3 .5*CR t IT ) S 251 

IT 3 1 S 252 

OETH 3 . 2*DE 1 H S 253 

NHH 3 1 S 254 

GO TO 3 S 255 

33 DR 1 1 ) 3 .2*DR (IT) S 256 

0ETH 3 .1*CETH S 257 

IT»1 S 258 

GO TO 3 S 2 59 

54 E4(IT) 3 ((RIC(4,NT1)/RW4C)-1.)*100. S 260 

IF (ABS(E4(IT)).GT.EM) GC TO 36 S 261 

NHC 3 1 S 262 

IC0UNT=IC0UNT*1 S 263 

Rw4C 3 RKBN S 264 

GO TO 42 S 265 
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C S 266 

C If BODY OThER THAN ft SPHERE IS CJNSIOEREO, APPROPRIATE CROSS SECTIO S 267 

C NAL BODY RADIUS AT JUNCTURE SHOULC BE USED. S 266 

C S 269 

35 CONTINUE S 270 

c S 271 

U THE FOLLCw I NG STATEMENTS COMPUTE CCNE CECMETRY. S 272 

C S 273 

ZC*ZIC(4»NT 1 ) - 22 S 274 

RT=ZC*TAN ( THC ) S 275 

Kw4C=Rl+RT S 276 

E4( I T) = ( (RIC(4»NT 1) /R W4C 1-1.01*100. S 277 

If i U.EU.l) EM* .001 S 278 

If ( ABS(E4( IT) ).L£.EM) GO TU 42 S 279 

36 CONTINUE S 280 

C S 281 

C THE FOLLOW I NG STATEMENTS COMPUTE CORRECTION TO RADIUS Of CURVATURE S 282 

C S 283 

If IIT.GE.2) GC 10 40 S 284 

IF IU.EG.il GC TO 37 S 285 

IF IE4I ITJ.GT.O.O) GO TO 39 S 286 

GO TO 38 S 287 

37 IF (E4IIT).GT.O..ANO.OR(IT).GT.O. I GO TO 39 S 288 

If (E4IIT).LT.0..AND.DRIIT).LT.0.I GO TC 39 S 289 

38 l T*l T+ l S 290 

DR(IT)*0R(IT-1)*2. S 251 

GO TO 3 S 292 

39 IT*IT+1 S 293 

OR 1 1 T)*DR( IT-ll-O. 50*ORI IT-1 ) S 294 

GO TO 3 S 295 

40 I T = I T+ 1 S 256 

10*10+1 S 297 

IF ( U.EU.l. AND. ID.GE.18.ANU.E41IT-1 ).LT..02) GO TO 42 S 298 

IF ( ID.GE .8 .AND . RKSX ( NO) .GT . 3000. > GO TO 41 S 299 

IF IRKSXINCI.GT. 10000.) GC TO 41 S 300 

IF <E4<IT-1J. EC. £4(iT-2>. ANO.NP.EQ.il GO TO 49 S 301 

If I IT. EG. 5 I £411) *E4( 7 1 S 302 

If IIT.EG.S) E4I2)*E4I8) S 303 

IF (IT.EQ.5) OR (2)*0RI8) S 304 

IF tIT.EQ.9l OR t 1 )=0RI 7 1 S 305 

IF t IT.EQ.5) I T*3 S 306 

OR! I T) *DR (IT-11— E441T-DMOR41 T— 1 l-CR( IT-2))/ tE4l IT-1 ) -E4 4 I T-2 ) ) S 3C7 

GO TO 3 S 308 

41 IF (RKSX(NC).GT. 10.) OR 1 1 ) *. 1*RKSX t NC ) S 309 

IF (RKSXtNCI.LT.10.) DR41)*-1. S 310 

OETH*-DTHC S 311 

RKSX4NC)*-RKSX(NC) S 312 

12*1 S 313 

IT*1 S 314 

GO TC 3 S 315 

42 CONTINUE S 316 

C S 317 

C THE FOLLOWING STATEMENTS CALL ROUT I RE TC WRITE CUNVERGEC RESULTS S 318 

C AND INITIALIZE CATA TO GO TO NEW POSITION ON SHOCK. S 319 

C S 320 

CALL OUTPUTS S 321 

If (RS(NO).GT.REND) GO TO 54 S 322 

If (THC. EG. 0.. ANC.IOPT.EQ.il GO TO 44 S 323 

IF (RHS.GT.O.) GC TC 44 S 324 

IF 4U.EQ.1) GO TO 46 S 325 

IF t THETASt NO) .LT. THSC) GC TO 43 S 326 

GO TO 44 S 327 

43 CALL FTLUP 4 THSC .RHS , 1 , NC, THETAS, RHC S ) S 328 

CALC FTLUP (THSC ,GSN, 1,NC» THETAS, US) S 329 

CALL FTLUP (THSC ,PSN,1, NC . THETAS , PS ) S 330 

44 IF (RKSX(ND).GT.IO.) GO TC 46 S 33i 


38 



APPENDIX B - Continued 


45 

46 


47 


46 

49 


50 


51 


52 


53 


IF ( ICCUNT.EC.4./ND. DETP.LT . . 5*CETF B ) DETH= .5*CETH6 
IE ( ICGUNT.GE.l.ANC. ICCUNT.LT.4) GO TO 45 
DKCli-.i 
GO TO 46 
DR ( 1 ) = DR( IT ) 

IF (RKSX(NDJ # GT.10.) DR (1) =0. 1 *RK SX ( NO ) 

DT HC*DET Ft 

IF l RKSXl ND) .LT.O. . AND. R K 3X { NU ) .GT 10 . ) CR (1 ) *-.15 
IF (RKSXIND) .LT.-10. ) CR l 1 ) = . 1*RK S X l ND ) 

IF (RKSX(NC ) .LT.1G0C. ) GC TC 47 

IF (IZ.EU.l) GC TO 47 

17*1 

DETH=-DETH 

RK$XIND)*~RK$X(ND) 

DR ( 1 ) *“DR ( IT ) 

IF ( IZ.EU.l. AND. A6S( THETAS (ND)-THSC) .LT..001) GO TO 49 
IF ( IZ.EU.l .ANO.(TFETAS(NC)-TH$C) .GT..001) GQ TO 49 
IF C IZ. EC. 1. AND. RKSXIND) .GT.O.) 0 RID = i*RKSX (ND ) 

IF ( IZ.EU* 1 -ANC.CETF.lt • * 001 ) DET H=- . i*UE TC 
IF (IZ. EC. 1. AND. RKSX (ND) .LT. -5000.) GO TO 48 
I T= 1 
NC*NG 
GO TO 2 

THS1=THETA$(N0) 

CONT INUE 
NC*ND 
N6*NC 
ZN=ZS(NC) 

RN=RS(NC) 

YC*Y 1 1 4t NT 1 ) 

RHOS(NB)*RhS 
US (N B) = USN 
PS(NB) =PSN 

IF (TFISI.CT .0. ) 1 F S* T HS 1 
IF ( THS1.EC.0.) lhS=TNSC 
XEFF-RN/SIN (THS) 

SI G* THS-THC 
YCC=XEFF*TAN(S IG ) 

YTM=YC-YGC 
ISUM^O 
I =NC + 1 
NZ*2 + 1 SUM 

IF (THS1.GT.0. .ANC.NZ.GT .5 ) T H$ = .5* ( THSC+THSi ) 

R S 1 = NZ *( . 15 ) 

RS( U =RN+RS1 

ZSl i I=ZN*(RS1/TAN( THS) ) 

SFS ( I j«0.5*RSm**2 

SIGMA=THETAIN-TFS 

OS Y* SFS (I )/M 1 

XS4 I ) * XEFF* (RSI/ SI N( THS) ) 

DO 5 i J= 1 1 NT 1 

SY( 4,J)=SFSCl)-J*DSY 

Y I (4, J )= ( 1.0/ (PFCS (NB ) *US ( N8 I *R S( I ) I ) * ( SF S ( I ) - SY ( 4, J ) ) 

CONTINUE 

Yl (4 f NTl)=Y I (4fNTl)+YTF 

YL = Y 1(4, Nil )/NTl 

DO 52 J = 1 tNT 1 

Yl ( 4, J)*J*YL 

CONTINUE 

DU 53 J* 1 *NT 1 

K 1 (4, J)*RS( I)-CGS(THS)*Y I(4 f J) 

RP*KS(I)-RI (4, J) 

ZI(4,J)=ZS< I)*(RP/TAN(SIGMA) ) 

PI 14 » J) =PS( NB) +( 1.0/ (RHOS(NB) *IRS (ll**2)*XS(I )**2))*(lSFSl I )**2)-S 
1 Y ( 4» J ) **2 ) 

CONTINUE 


S 332 
S 333 
S 334 
S 335 
S 336 
S 337 
S 338 
S 33S 
S 340 
S 341 
S 342 
S 343 
S 344 
S 345 
S 346 
S 347 
S 346 
S 349 
S 350 
S 3 5 1 
S 352 
S 353 
S 354 
S 355 
S 356 
S 357 
S 358 
S 3 59 
S 360 
S 361 
S 362 
S 37 3 
S 364 
S 365 
S 366 
S 367 
S 36b 
S 369 
S 3 70 
S 371 
S 372 
S 373 
S 374 
S 375 
S 3 76 
S 377 
S 378 
S 379 
S 380 
S 381 
S 362 
S 363 
S 384 
S 3 6 5 
S 386 
S 387 
S 388 
S 389 
S 390 
S 391 
S 392 
S 393 
S 394 
S 395 
S 396 
S 3S7 
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KGC ( 4 ) *RS ( I ) S 398 

ZCC(4)=ZSU) S 399 

RIC(4,NTl)*RI(4tMl) S 400 

YIC(4fNTll*YI(4*MlJ S 40X 

Z I C ( 4fNTli*ZI ( 4 *NT li S 402 

CALL OUTPUTS $ 403 

IF (RSUi.GT.RENO) GO TO 54 $ 404 

I SUM= 1SUM+X S 405 

NC*I S 406 

GO TO 50 S 407 

54 RETURN S 408 

C S 409 

55 FORMAT (//2X»3XH BODY CONE JUNCTURE IS LOCATED ) S 4X0 

end s 4XX- 


c 

c 

c 


X 


c 

2 

3 


SUBROUTINE OUTPUTS 

THIS SUBROUTINE WRITES T PE SOLUTIONS FOR SUPERSONIC FLOW FIELD 

COMMON /BLK 2/ IT »OELX, I BEGf EMtNCt NT X t NX»REND 9 ND«PST«uR( 8 ) «F (200) t 1 
XP 

COMMON / BLK4/ XS ( 300) » RKSX( 300) f Z Si 3 00) »RS ( 300) t SFS ( 300) t TP ET AS( 30 
X 0) 9 US( 300) tPS(300) fRPiOS (300) t MS (300 ) * VS ( 300 ) 

COMMON /BLK 5/ R I ( 7 » X 5) » L I ( It X 5 ) » Y I ( 7 . X 5 ) f R I C( 7 , X5 ) » Z IC ( 7 1 X5 ) t Y IC (7 
X tX5) t RW4C t ZSTtE4(8)vTHB(50) 

COMMON /BLK6/ SFC( 7) 9 RO 4 7) , ZO i /) , PO ( 7) 9 U0( 7 ) , THE T AO ( 7 ) f RHOSC ( 7 ) * PS 
X 0 ( 7 ) ,RK0(7> ,DELSY( 7) f RKOC47 ) 9 ROC ( 7 ) » ZOC ( 7 ) » VO ( 7) 

COMMON /BLK 7/ SY ( 7 t X 5) 9 P I ( 7 , X 5) , RHI I 7 , X5 ) ,U 1 ( 7 9 X5 ) » MS 1 4 7 9 X5 ) » V I ( 7 * 
X X 5) 

IF (IC.EU.2) GC TO X 

IB-NX-X 

ib*ia*x 

IC S 2 

WRITE (6 9 2) ZIC(4,NTX ) 9RIC(4»NTX) , Y I C ( 4 , NT X ) , ZCC ( 4) , ROC (4 ) 9 NO 

WRITE (6 t 3) PS (ND) 9 RHGS ( NO ) t US ( NO ) 9 P I ( 4t NTX) , RHI ( 4 f NTX ) # UI ( 4. NTX ) 9 
lMSI(4 f NTX)*THB(ie) 

RETURN 

FORMAT (//2X93HZW=9EIX*4 f 2X»3HRW*,EXX.4.2X»3HYW*»EXi*4,2X93HZS*.EX 
X X* 4 #2X , 3PiRS* 9EXX«4tI3) 

FORMAT ( 2X 9 3HP S* 9 E X X . 4 f 2X , 5HRH0S** # E X X . 4 »2X f 3HUS= 9 EX 1.4 9 2X 9 3hPW= * E X 
XX.492X95HRHCW**EXX.4,2X93hUW=»EXi .492X t 3HMWs, ElX.4 f /2X f AHThfl^ # EXX. 
24) 

END 


T X 
T 2 
T 3 
T 4 
T 5 
T 6 
T 7 
T 8 
T 9 
T XO 
T XX 
T X2 
T 13 
T 14 
T X 5 
T X 6 
T 17 
T X8 
T X9 
T 20 
T 21 
T 22 
T 23 
T 24 
T 25 
T 26 
T 27 
T 28 
T 29- 


The sample output is for a paraboloid at = 10 and y =1.4. 


STAG P T SHK RAO CUR - 1.42166*00 STAG PT STAND OFF 01 ST - 1.47436-01 

flOUr STAG PT PRESS* 9.2 2 9dE-Q 1 


ARC LENGTH 

SHOCK AXIAL LCC 

SHOCK RADIAL LOC 

BODY AXIAL LOC 

BODY RADIAL LOC 

1 * 50 77E- 01 

7.9416E-03 

1 • 5049E—01 

1.56 1 6E— 0 1 

1. 3463E-01 

3.0155E-01 

3. 1465E-02 

2. 99356-01 

1.83236-01 

2.6/56E-01 

4 . 5232E-01 

6. 9379E-Q2 

4.4522E-01 

2 . 260 66— 01 

3.S670E-01 

6.0309E-01 

1.201 36-01 

5.671 5E-01 

2. 8 220E-01 

5.21CBE-Q1 

7.5JB7E-01 

1.82016-01 

7.2460E-01 

3.5245E-01 

6. 4022E-01 

9.04646-01 

2. 53356-01 

8. 5 739E-01 

4.31 746-01 

7.540 76-01 


BUOY NORMAL LOC 
4699E-01 
1.55046-01 
1.64C2E-01 
1.75S5E-01 
1.9019E-01 
2.0616E-01 


SHOCK RADIUS OF CURV 
1.45 136*0 C 
1 • 53956 *0 C 
1.6B656*0 C 
1* 69226*0 C 
2. 15676*00 
2.47996*00 
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XSP- ..b23<IE-01 


7.6636E-01 

5. 7C50E *00 

2.90706-01 

1.58706-01 

7.63716-01 

9.25766-02 

7.72876-01 

5. 71766+00 

2.85506-01 

1.4926E-01 

7.40566-01 

8.5964E-02 

7.78706-01 

5.72626+00 

2.8032E-01 

1 . 3960E— 0 l 

7.17716-01 

7.9351E- 02 

7.83866-01 

5.73C3E+00 

2.75166-01 

1.2968E-01 

6.95146-01 

7.2 739E-02 

7.88346-01 

5.73 156 +00 

2.6999E-01 

1 .19506-01 

6.7284E-01 

6.6126E-02 

7.9215 E— 01 

5. 72846+00 

2.64 79E-01 

1.09026-01 

6.50816-01 

5.95136-02 

7 .9 52d t-Ol 

5- 7215E+C0 

2.59526-01 

9.82316-02 

6 • 2903E-01 

5.29016-02 

7.97746-01 

5*7 1QSE +00 

2.54156-01 

8.7056E-02 

6.07516-Cl 

4.6288E-02 

7.99536-01 

5.6S66E+0Q 

2.4d636-01 

7.55956-02 

5.862 5 E-0 1 

3. 9 6 76E- 02 

8 .00636—0 1 

5.67866+00 

2. 42926-01 

6.37056-02 

5 . 6525E-CI 

3.3063E-02 

8.0107E-01 

5- 65 716+CC 

2.36936-01 

?.141 16-02 

5.44526- Cl 

2.64506-02 

8.00826— 01 

5 • 63216+00 

2.30606-01 

3.871 66- C 2 

5.24076-01 

1 .98386-02 

7. 99916-01 

5. oC 356 + 00 

2.2361E-01 

2.5664 t-Oc 

5.03926-01 

1.3225E-02 

7.93J1E-01 

5*57156+00 

2. 16466-01 

1.24256-02 

4.8410E-01 

6.6126 E-03 

7.96056-01 

5. 53616 +0 C 

2.0847E-01 

0 . 

4.64636-01 

4.4409E-16 


MACH NO ML 

SHOCK PRESSURE, 

TAN SHK VELOCITY 

SCR SHK VELOCITY 

ShCCK OENSITY, 

SHOCK MACH NJMd 

80 OY AN6L! 

1.7362E-01 

8.22906-01 

1.05126-01 

1. 74126-01 

5.71126+00 

4.5304E-01 

i, 43686+00 

3.21866-01 

7. 57066-01 

2.05026-01 

1. 71646-01 

5. 70236+00 

6.0457E-01 

1.30946+00 

4. 6463c - 01 

7. 591 dE -01 

2.55626-01 

1.67936-01 

5 . 6d£3E+00 

7.tt7C2E-0i 

1.1931E+0C 

5.96586-01 

7. 1463E-01 

3. 15476-01 

1. t34fa6-0l 

5.67006+00 

9. 7454E-01 

1.09046+00 

7. 15/66-01 

6.6d09E-0! 

4.43676-01 

1.58666-01 

5.64846+00 

1. 1579E+00 

1.0013E+00 

8.21806-01 

6.22806-01 

5.01216-01 

1. 53656-01 

5. 62446 + 00 

1.53166+00 

9.2470E-0 1 


XSP- 1.5077E-01 


FIELD PRESSURE 

FIELU DENSITY 

FIELD VELOCITY (TANS) 

FIELD VELCC ITY(NGRK) 

FIEtO MACH NO 

FIELD STRM FUNC 

8.3366E-01 

5.761 7E+00 

1.0 1 5oE-0 1 

1.63666-01 

4. 28326-C1 

l. 06126-02 

8.43666-01 

5. 8C616+CC 

5.833 3E-02 

1.5342E-01 

4.04096-01 

9. 8537E-03 

8. 52d9 E— Oi 

5.85056+00 

9.56466-02 

l. 42756-01 

3.80436-01 

9.095 7E-03 

6.61356-01 

5. dd90t +0C 

9.34516-02 

1.31576-01 

3.5 7356- Cl 

8.337 7E -03 

8. 69046-01 

5.92366+00 

9. 16SoE-02 

1.20956-01 

3.34906-01 

7.5 757E-03 

a. 7597E-01 

5.5543E+00 

9.03166-02 

1.09726-01 

3. 13136-01 

6.82186-03 

8.82 126-01 

3.58 126 +00 

8.92376-02 

9.8265E-C2 

2.92126-01 

6.06386-03 

8.8 7526-01 

6.0C43E+0C 

8.83766-02 

8.65836-02 

2* 7 IS 76 -Cl 

5. 3058E- 03 

8.92146-01 

6.0236E +00 

6. 76 386-02 

7. 46606-02 

2.52836-01 

4.5478E-03 

8.9600E-01 

6.C392E+00 

d. 64176-02 

6.2490E-02 

2 • 34896- C 1 

3. 7899E-03 

8.990a 6-0 1 

6.05116+00 

8.605 76-02 

5.00756-02 

2.18386-01 

3.0319E-03 

9.0141E-01 

6.05926+00 

8.50496-02 

3. 74286-02 

2.03ol£-0l 

2.27396-03 

9.029©E-0i 

6. 06 3oE +00 

8.36386-02 

2.46166-02 

1.9095E-01 

1.51596-03 

5.03756-01 

6. C6446+00 

8.1 73 7E-02 

1. io20E-0^ 

1.80816-01 

7.5757E-04 

9. 03766-01 

6.06 146+00 

7.9323E-02 

0. 

1. 7362E-01 

5.55116-17 


FIELD NORMAL LQC 

FIELD AXIAL LOC 

FIELD RADIAL LOC 

8.52316-03 

1.64fciE-02 

1.49596-01 

i. 73186-02 

2.52136-02 

1.48b 76-01 

2. 6 368 E- 02 

3.42146-02 

1.4772E-01 

3.50486-02 

4.34426-02 

1 • 4b7 46- 0 1 

4.51316-02 

5.28726-02 

l .437 56-01 

5.47896-02 

6.247 86-02 

1. 44736-01 

6.46016-02 

7.22356-02 

1.4370E-01 

7.45506-02 

8.2129E-02 

1.4265E-01 

8. 46276-02 

5.21496-02 

1 .41596-01 

5.48346-02 

1.02306-01 

1. 40p 26- Cl 

1.0519E-01 

1. 12&0E-01 

l . 394 36— C 1 

1.15726-01 

1.23076-01 

1 * d 83 36- C 1 

1.26476-01 

1.33766-01 

A. 37206-01 

1.3753E— 0 1 

1.44766-01 

l . 3603E-0 1 

1 .48996-01 

1.5616E-01 

l .34836- G 1 


XSP« 3.01556-C1 
8.0630E-01 

5.73896+00 

2. 00726-01 

1.61S5E-01 

5.81526-01 

4.18866-02 

8.l4d06-01 

5. /713E+OC 

1.96476-01 

1. 52076-01 

5. 5883fc— 01 

3.6899E-02 

8.225/ E-Oi 

5. 7557E+0C 

1.92426-01 

1. 41986-01 

5.3,6656-01 

3.59026-02 

8.29O2E-01 

5.82426+00 

1.88576-01 

1. 3 1666-01 

3. 15006-01 

3.291 06-02 

8.3594 E—01 

5 • 84 4 86 + 0 C 

1.848 7E- 01 

1.210SE-01 

4.93886-01 

2.99186-02 

8.41536-01 

5. 86166+00 

1 .d L316-01 

1. 1025E-01 

4. 73326-01 

2.69266-02 

8 .48396—01 

5 . a 7 456 +0 C 

1. 7765E-01 

9.91226-02 

4.5335E-01 

2.3935E-02 

8.505^6-01 

5.86386+00 

1.74446-01 

8 . 7bdSE- C c 

4.34006-01 

2.09486-02 

d. 53936-01 

5.88536+00 

1. 7103E-01 

7.59276-02 

4.15336-01 

1.795 IE- 02 

8.5660E-01 

5. d9 1 1 E +00 

1.6755 E-01 

6.38226-02 

3.973SE-01 

U4959E-02 

8.5655E-Q1 

5. 86526+00 

1.6393E-01 

5.1d6oE-02 

3 . 6026E-G 1 

1. 19676-02 

8. 5977 E— 01 

5.8837E+00 

1.60C7E-01 

3. 8570E-C2 

3* 64036- C 1 

8.97556-03 

8.602b E-01 

3. 6 7466 +C C 

1.55866-01 

2.54d7E-02 

3.48 606— Cl 

5. 98366-03 

6.6002 E-Oi 

5.86196+00 

1. 51196-01 

1.2297E-0* 

3*34 706-CL 

2.99186-03 

8.5905E-01 

5. 84576+OC 

1. 45996-01 

0 . 

3.2186E-01 

2.2204E-16 
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APPENDIX B 
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FIELD NOR MAI LUC 

FIELD AXIAL LOC 

FIELD RADIAL LQC 


8. 63666-03 

3. 993 8E -02 

2 . 9 75 86- Cl 


1.74606-02 

4 • 857 4E -02 

2.95776-01 


2.64816-02 

5.7403E-02 

2.9392E-01 


3.5705E-02 

6.643 IE- 02 

2.92036-01 


4.5 1406-02 

7.56666-02 

2.90106-01 


5.47976-02 

8.5U8E-02 

2.88126-01 


6. 46856-02 

9.4796E-02 

2.86096-01 


7.48206-02 

1.04726-01 

2.840 IE-01 


8.52196-02 

i. 64896-01 

2 .8 Id 8 6-0 1 


9.59046-02 

1. 253 56-0 l 

2. 79696-01 


1.0690E-01 

1.3612E-01 

2.77436-01 


1.18266-01 

1.4723E-01 

2.75116-01 


1 • 3002 6-01 

1.58746-01 

2.72706-01 


1.4225E-01 

1.70716-01 

2.7019E-01 


1. 55046-01 

1.83236-01 

2.67576-01 



FIELD NORMAL LOC 

FIELD AXIAL LOC 

FIELD RADIAL LUC 

6. 9 lbOE— 03 

7. /d 966 -02 

4 .425 9 E-01 

1.O024E-02 

b . 65966- 0 2 

4.3989E-C1 

2.73396-02 

9.54946-02 

4.3714 E— 01 

3.68776-02 

1.046 16-01 

4. 343 2E - 01 

4. 6653 E— 02 

1. 15956-01 

4.314 2E- 0 1 

5.6694E-02 

1. 2 3546 ~ 01 

4. 2d45 E- C 1 

6.70L3E-02 

i • 3339E-0 1 

4 . 2540E- 0 l 

7.76406-02 

1. 4354E-01 

4 • 222o£- 01 

8.b603 6-02 

1.5402E-01 

4. 190 IE- 0 1 

9. 9937E— 0^ 

1.64846-01 

4.1566E-01 

1.1 lobE-0 1 

1. 7606E-01 

4.12196-01 

1.23396-01 

1.87736-01 

4.085 7E-01 

1 .36636-0 1 

1. 9989E-0 1 

4.046 IE-01 

1. 49976-01 

2 . 1 26 46 - 0 1 

4.00866-01 

1.6402 E-01 

2.26C6E-01 

3.96706-01 


XSP- 6.0309E-0L 
7. 19666-01 

5. 66 446*00 

3.o975 E-01 

1.54716-0 1 

9.5035E-01 

1. 6094E-0 1 

7.24086-01 

5 • 65 46E *00 

3.63656-01 

1. 45746-01 

9.25726-01 

1. 49446-01 

7.27896-01 

3.6407E+0C 

3.57866-01 

1 . 3o5 3 E-0 1 

9.01126-01 

1.37956-01 

7.31076-01 

5 . 622 BE +00 

3.51 73E-01 

1.27056-01 

8.76546-01 

1.2645E-01 

7.33656-01 

5.60106+00 

3.45456-01 

1.17296-01 

8.51936-01 

1.14966-01 

7.3561E-01 

5.57 53E+0C 

3.59006-01 

1.07226-01 

8.27276-01 

1.03466-01 

7. 36956-01 

5.54586+00 

3.5233E-01 

9.68046-02 

6.02526-01 

9.19666-02 

7.37686-01 

5.5127E+CC 

3.25416-01 

8.60216-02 

7.77656-01 

8. 04706-02 

7.37/96-01 

5.47586+00 

3.18196-01 

7 • 46 4CE-02 

/. 52616-01 

6.89756-02 

7.37296-01 

5 • 43 546 +C C 

3. 106 IE- 01 

6.32326-02 

7.27366-01 

5.74796-02 

7.3617 E-Oi 

5.89146+00 

5.02596-01 

5.1173E-C2 

7.C191E-C1 

4.59836-02 

7.34446-01 

5.54406+00 

2 ■ 94GoE - 01 

3.86556-02 

6.76156-01 

3.44876-02 

7.3209 E-01 

5.29516+00 

2.8489 E- 01 

2.571C6-02 

6.50056-01 

2.2992E-02 

7.29136-01 

5.2389E+CC 

2.7 49 b6~ 01 

1 .24946-02 

6.23556-01 

1.14966-02 

7.2 555 E-01 

5. 18136+00 

2 • 64 15 E- 0 1 

0. 

5. 96586-01 

a. 6818 6- 16 


FIELD NOR HAL LOC 

FIELD AXIAL LOC 

FIELD RACIAL LOC 

9.29866-03 

1.2874E-01 

5 ■ 836 66- C 1 

1. O815E-02 

i. 3 7566-01 

5 • dOO 96-01 

2.85596-02 

1 ■ 4660E-0 1 

5. 76426-01 

3.85836-02 

1.5589E-01 

5.72566-01 

4.8876 E-02 

I.6543E-01 

5.68806-01 

5.9483E-02 

i . 7 52 6t- 01 

5.64826-01 

7.04266-02 

1.854DE-01 

5.607 16-01 

8.17436-02 

X . 958 9E-0 1 

5 .56466- 0 1 

9.54736-02 

2.0o7t£-01 

5 .52 05 E- 0 1 

1.0566 E-01 

2. 18066-01 

5.47486-01 

1.18 37E-Q1 

2. 2 983E-0 1 

5.427 16-01 

1.3lb6 E-01 

2. 42156-01 

5.37726-01 

1.4561E-01 

2.550BE-01 

5.3248E-G1 

1.60336-01 

2.6d 72E-01 

D . 269 56- C 1 

1.7 5956—0 i 

2.8320E-01 

5.21096-01 


XSP= 7.53876-01 
6 . 7 1 1 7 E-0 l 

5 • b2 C 7E *00 

4.3761 E-01 

1.50396-01 

1.1317E+00 

2.45066-01 

6. 7 3 70E-01 

5. 5fc 876 + 00 

4.31526-01 

1 .41 88 E-0 1 

1.10536+00 

2.2755E-01 

6. 756 7 E— 0 1 

3 • 3525E *00 

4.2481 E-01 

1.35136-01 

1.07 866 + CO 

2.10056-01 

6.77096-01 

5. 5 1 2 it *C C 

4.1a 056-01 

1.24116-01 

1.0516E+00 

1.9254E-01 

6.7 795 E-Oi 

5 . 46 8 IE +00 

4.11026-01 

1.14786-01 

1. 02436+00 

1.7 504E-0 1 

6. 7826E-01 

5.41996+00 

4. 036 7E-0 i 

1.05135-01 

9.9658E-C1 

1.5754E-C1 

6.7801 E-01 

5.56796 +00 

5.9 59 7E- 0 1 

9.51186-02 

9.68416-01 

1 .4003 E-01 

6.77^06-01 

5. 5 1 20E +00 

3.87856-01 

8. <*7156-02 

9. 39 7 IE- 01 

1.22536-01 

6.75846-01 

5. 25 256 +0 C 

3. 792 7E- 01 

7.38866-02 

9.10416-01 

1.05026-01 

6. 7393E-Q1 

5. 18946+00 

3.7015E-01 

6.25946-02 

8.80416-01 

8. 75206-02 

6.7 1466-0 i 

5. 12266+00 

3.60396-01 

3 • 08C8E-02 

8.4962E-01 

7.00166-02 

6.6843 E-01 

5.05246*00 

3.49896-01 

3.8 50 7E- 02 

8. 17916-01 

5.25126-02 

6.64856-01 

4 . 9 7 ddt +00 

3. 38506-01 

2.3/0/6-02 

7.85146-01 

3.50066-02 

6.607cE-01 

4. 9018E+00 

i • 2bObE- 0 1 

1.25466-02 

7.511 6£ -Cl 

1.75046-02 

b. 56026-01 

4. 02 16E+00 

3.1239E-01 

0 . 

7. 15766-01 

1.77646-15 


[ 
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APPENDIX B - Concluded 


FIELD NORMAL LUC 

FIELD AX I 3 L LUC 

FIELD RACIAL LOC 

4.75796-03 

1.90756-01 

7.202 76-01 

1.97686-02 

1.99726-01 

7. 15d it-Gl 

3.0C556-02 

2.08946-01 

7.11266-01 

4.06446-02 

2. 18436-01 

7.0657E-G1 

5.1 5666—02 

2.2822E-01 

7.017 2E-01 

6.28546-02 

2. 38346-01 

6 • 9b7 16-0 l 

7.4546E-02 

8.48816-01 

6 .91326-01 

8.66876-02 

2. 547 06 -Cl 

O. 86146-01 

9.93296-02 

2.71026-01 

6.80536-01 

1.12536-01 

2 . 82866-01 

6. 74676-01 

1.26376-01 

2.95266-01 

6 .68536-01 

1.40436-01 

3.08316-01 

6.6207E-01 

1.56326-01 

3.22106-01 

6.5524E-01 

1.7268E-01 

3. 3o 766- 0 1 

6.4798E-01 

1. 90196-0i 

3.5245E-01 

6.40226-01 


XSP- 9.04646-01 
6.242 76-01 

5.57676*00 

4.95056-01 

1 .46026-01 

1.30386*00 

3.43066-01 

b. 25246-01 

5.52456*00 

4.88616-01 

1,379/6-01 

1.27556*00 

3.1856E-01 

6.25726-01 

5.4t6l£*0C 

4. 81856-01 

1.29656-01 

1 .24676*00 

2.94056-01 

6.25696-01 

5. 40 /4E *00 

4.74756-01 

1.21066-01 

1.21 736*00 

2.69556-01 

6.25166-01 

5. 34 26t *C 0 

4.67256-01 

1.12186-01 

1.18726*00 

2.45046-01 

6.24136-01 

5.27376+00 

4.59326-01 

1 . 0292E-0 1 

1.15646*00 

2.20546-01 

6.22616-01 

5. 20C9E+00 

4.50896-01 

9.33146-02 

1.12476*00 

1.96036-01 

6.20586-01 

5.12406+00 

4.41916-01 

8.32956-02 

1.09216*00 

1.71536-01 

6.16056-01 

5. 04346*00 

4.32276-01 

7.28276-02 

,1.05836*00 

1.47036-01 

6. 15036-01 

4. 958 9E *00 

4.21906-01 

6. 1866E-02 

1.02336*00 

1.22526-01 

6.11506-01 

4. o707£*C C 

4.10666-01 

5.03716-02 

9.66876-01 

9.80176-02 

6.07476-01 

4.77906*00 

3.98416-01 

3. 830 7E-02 

9.48796-01 

7.35136-02 

6.02956-01 

4. 6fc 3aE *0 0 

3 • 845 b E— 0 1 

2. 56 756-02 

9.08816-01 

4.90086-02 

5.4 792 E— 0 1 

4. 66496+OC 

3. 70086-01 

1 .25876-02 

8 .60 bl E-Ql 

2.45046-02 

5.9239E-01 

4.48276*00 

3. 53486-01 

0. 

8.21 606- Q 1 

1.77646-15 


FIELD NORMAL LUC 

FIELO AXIAL LUC 

FIELD A AD 1 Al 

1 . 0 c74E-02 

2. 62246 -0 L 

o. 52256*01 

2.08396-02 

7 138E ~ C 1 

8 . * d9t»6- Cl 

3.17216-02 

2.80806—01 

8.4i30e-Ci 

4.29546-02 

2. 90526- Cl 

8.3587E-01 

5.457j 6-0*: 

3.00576-01 

8.3004E-C1 

6 . 662 1 6- 02 

3. 11006-01 

8.24006-01 

7.91*46-02 

3.2ldJE-0i 

8.17736-01 

9.21986-0*: 

3.33136-01 

a. in dE-oi 

1.05856-01 

3. 44946-01 

8.0*346-01 

1 . 20 L 76-01 

3.57346-01 

7.97loE-0i 

1.35276-01 

3. 70*06-01 

7. 89606-01 

1.5125E-0I 

3 • 842 3c- 01 

7.81596-01 

1.6626E-01 

3. S8946-01 

/. /306E-01 

1.86486-01 

4. 14716-01 

/ .6393E- C i 

2.06166-01 

4. 31746-01 

7.340 76-01 


ZW = 5.0366E-01 = 8.4*076- 01 Yrf* 2.20556-Ui I 5 = 3.18346-01 R S* 9.6364t~Oi *3 

PS* 5.87206-01 RHUS* 3.oC326*Q0 U$* 5.4215E-01 PW* 3.42986-01 RhUw* 4.2124E*00 Uw = 3.843C6-QI Mw= 9. 04636-01 

THB* 8.69 7 5E-0 1 


Zw* 5.8532E-01 KM* 9.3583E-01 YW = 2.36586-01 ZS = 3.92756-01 RS* l. 0732E+Q0 5G 

PS* 5.509/6-01 KHOS* 5.57908*00 U$* 5.80876-31 PM* 4.936/6-01 RHU«* 3.93546*00 Uw= 4.146GE-01 M* = 9.89346-01 

THB* 8.18536-01 


ZW* 6. 7 6976—0 1 Rw* 1.02916*00 Y** 2.54126-01 l S* 4.7687E-01 R5* 1.1858E+Q0 57 

PS* 5.15486*01 khGS* 5.55226*00 US* b.io446-0l P** *.5C47E~0i AHQ** 3.68916*00 LW= 4.40806-01 MW* 1.06556*00 

THB* 7.7 10 5E-01 


ZW* 7.8 124E-0 1 R** 1. 12 5 9£*C0 YM* 2. 72616-01 ZS* 5.73116-01 AS* I.i035E*00 64 

PS* 4.8104E-01 RHUS- 5.5228E*00 US* 6.4909E-01 PM* 4.1Q71E-0L RHCJ** 3.45086*00 CM* 4.65696-01 MW* 1. i40aE*C0 
THB* 7.26256-01 


Zw* 8.98626-01 8** 1.22576*00 Y«* 2.95C66-01 ZS= 6. 82006-01 AS- 1.42616*00 7 1 

PS* 4.4799E-01 RHUS* 5.49C7fc*Q0 uS= 6.78966-01 P#» 3.7355E-01 RhUw* 2.2248E*00 UW* 4.89016-01 Mw* 1.21436*00 

THB* 6.84336-01 


ZW* 1 . 029aE*00 RW* 1.32856*00 YM* 3.1d53E-01 L S= 8. 04326-01 RS- l. 5534E*00 78 

PS* 4. U52E-01 RHUS* 5.45586 + 00 US* 7.0622E-01 Pit* 3.3936E-01 RHOm* 3.01116*00 Uw* 5.10956-01 Mm* 1.28636*00 
THB* 6.4526E-01 
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(a) Shock shape. 
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Ref. 7 

Present results 











(a) Shock shape. 



(b) Body pressure distribution. 

Figure 9. - Comparison of shock shape and pressure distribution on 
hemisphere-cylinder at Mqq = 7.7 and 7 = 1.4. d = 0.75. 



(a) Shock shape. 



(b) Body pressure distribution. 

Figure 10. - Comparison of shock shape and pressure distribution on spherically 
blunted 30° half-angle cone at - 8 and 7-1.4. r n = 1.0. 






